1. Introduction

We take a model from literature and look at its major characteristics. Then we ask three questions about it:

  • How does the model output respond to a single pulse perturbation?
  • How does it respond to repeated (periodic) perturbations?
  • How does the model dynamics impact on downstream connected elements?

The answer to the first question will tell us what the model’s dynamical repertoire is, both in parameter space and in state space. It is at the heart of the “encoding” problem.

The answer to the second question will tell us how the model’s dynamics change in response to time-dependent input from the environment. It touches on the problem of interaction between time scales in different cellular components.

The answer to the third question will tell us what the model’s dynamics do to other functional components, i.e. how it connects to other dynamical elements of a cell. It is at the heart of the “decoding” problem.

Overall this contributes to the broader question: how does a signalling network respond to time-dependent input and how is this response converted to reliable downstream functionality.

2. Biological Finding and Research Question

We examine a published model of calcium signalling to account for experimental observations of cytosolic calcium levels in response to agonist stimulation.

The experimental finding that started the process of modelling was the following measurement of cytosolic calcium in isolated hepatocytes under stimulation with ATP, UTP or ADP:

image number 1

This shows that the same intracellular signalling pathway may respond with different types of dynamics to different external stimuli. The pattern observed with ADP is referred to as regular spiking. There is a preferred frequency and a fairly constant amplitude of the spikes throughout the period of stimulation. A similar response is seen in hepatocytes when treated with constant levels of phenylephrine (see below). The pattern observed with ATP/UTP, in contrast, is referred to as bursting. It is characterised by a main large-amplitude spike followed by a number of secondary oscillations of small amplitude at elevated mean level of calcium and a spontaneous return to basal levels. This process self-repeats when the stimulus continues and stops when the agonist is washed out.

Previously published models of calcium signalling showed periodic spiking but were unable to account for the more complex bursting pattern. It is therefore crucial to have a model that can display both types of responses.

The model assumes agonist stimulation to lead to the dissociation of the G-alpha subunit of a receptor bound G-protein.This activated G-alpha subunit is then assumed to activate membrane bound protein phospholipase C ($PLC$) and lead to the release of calcium ions from intracellular calcium stores in the endoplasmic reticulum (ER).

Schematic of the induction of cytosolic calcium fluctuations when signalling is started via activation of a G-protein (red and green) bound membrane receptor (purple). Protein kinase C (PKC) is one of many downstream proteins whose activity depends on the presence of calcium ions.

We will first describe the model and use basic mathematics and MATLAB® to study its main features. Then we will consider perturbations to the model in the form of time-varying agonist concentrations. Finally we will look at the downstream effects of different type of calcium dynamics on calcium-dependent proteins, like PKC and other calcium-binding proteins (CBPs).

image number 2
3. The Data

A look at the experimental data with the corresponding Fourier spectrum.

image number 1
image number 2

The frequency $f$ the regular spiking is around $2.5 ~min^{-1}$. Some harmonics can be see at around $5$ and $7.5 ~min^{-1}$. This indicates that the signal is not pure sinusoidal. Bursting is slower with a peak at around $0.4 ~min^{-1}$. There are also some harmonics but the signal is too noisy to give clear peaks in the higher frequencies. It is worth noting that the number of secondary oscillations varies leading to variation in width of the bursts.

A brief tutorial on how to import the data to MATLAB® and calculate their Fourier spectrum is provided separately.

4. The Model

Kummer et al. (2000) developed a minimal model where simple spiking and bursting can be studied. Here are the equations:

$$ \begin{align} rate(G_\alpha) = & k_0 + k_1 \times G_\alpha \\ & - k_3 \times G_\alpha \times \frac{PLC}{k_4 + G_\alpha} \\ & - k_5 \times G_\alpha \times \frac{C_\alpha}{k_6 + G_\alpha} \end{align} $$
$$ rate(PLC) = k_7 \times G_\alpha - k_8 \times \frac{PLC}{k_9 + PLC} $$
$$ rate(C_a) = k_{10} \times G_\alpha - k_{11} \times \frac{C_a}{k_{12} + C_a} $$

$G_\alpha$ denotes the alpha subunit of the receptor-coupled G-protein; variable $PLC$ represents phospholipase C; and variable c represents cytosolic calcium. Other players (like IP3 and stored calcium in the ER) are taken into account only implicitly.

According to this model, there there are two “up” processes that contribute to the presence of activated $G_\alpha$ protein: a constant activation ($k_0$) which can be tuned to account for different levels of agonist in the environment; and an activation that depends on the $G_\alpha$ protein level ($k_1 \times G_\alpha$). And there are two “down” processes: a decay dependent on $G_\alpha$ protein and $PLC$ (controlled by $k_3$ and $k_4$); and a decay dependent on $G_\alpha$ protein and $C_a$ (controlled by $k_4$ and $k_6$).

The rate equations for $PLC$ and $C_a$ each have an “up” process dependent on the $G_\alpha$ protein level ($k_7 \times G_\alpha$ and $k_{10} \times G_\alpha$, respectively); and a “down” process depending on their own concentration in a Michaelis-Menten fashion controlled by constants $k_8$ / $k_9$ and $k_{11}$ / $k_{12}$ , respectively.

The model forms a reaction network composed of three nodes and eight (inter-)actions involving nodes.

The appendix gives numerical values for the model constants which are taken from the original publication.

Note that this is a simplified version of a more detailed model. The simplified version was shown to qualitatively reproduce dynamical features of the detailed model.

Mathematics: Rate Equations

Rate equations belong to the class of differential equations. In (bio-)chemical reaction kinetics the variables are often expressed as concentrations, i.e. number of particles (molecules or ions) per volume or area. Then the rate equation is an instruction how to calculate the rate of change of a given concentration over time. In in vitro kinetics it is often assumed that the concentration of substrates or products of a reaction changes as a function of only a single concentration. However, in biochemical networks any concentration typically depends on more than one concentration. This is also the case in our model of calcium signalling.

As a consequence, the rates of changes of all concentrations involved in a network need to be considered at the same time. We then do not deal with a set of independent equations but with a set of interconnected differential equations, in our case 3 coupled differential equations. In the form we wrote them they are called “ordinary” differential equations.

Note that the current model does not include any assumptions about e.g. cell morphology or about the fact that G-protein subunit alpha is membrane bound, whereas calcium ion concentration changes throughout the cytosol of a cell. This is a simplifying assumption and is made mostly for convenience. For example, if the spatial distribution of cytosolic calcium were to be described by the model, it would have to be of a more complicated type that describes concentration changes in different locations. The resulting equations, so-called “partial” differential equations, are mathematically more difficult to treat and demand more computational power to simulate.

5. Simulation of the Unperturbed Model
5.1. Regular spiking

If the agonist concentration is increased from 0 to a constant level (bottom) for small to medium values of $k_2$, regular spiking results (top). With further increase of agonist levels (increasing $k_1$), the frequency of the spiking increases while the amplitude stays relatively constant (amplitude coding).

image number 1
MATLAB®: How to create the agonist protocol graph

The bottom graph in the above figure shows the shanges in parameter $k_1$ over time. An array of values for $k_1$ to be plotted in this way can be set up using the function “repmat”:

% Base value of k_1
k_1 = 0.01;

% Produce an array with 1000 copies of k_1
k_1_array = repmat(k_1,1000,1)
% Change a subset of these values("Stimulation")
k_1_array(500:600) = 0.2;

% Plot graph
5.2. Regular bursting
image number 2

At higher values of $k_2$ (to account for modified receptor kinetics) regular bursting is observed. Within the bursting dynamics, frequency does not change much if $k_1$ is varied.

6. Analysis of Model Structure

The core model contains only a small number of terms. This make an analysis of the model feasible.

Note the following features:

  • The equation for $G_\alpha$ depends on variables $G_\alpha$, $PLC$, and $C_a$.
  • The equation for $PLC$ depends on variables $G_\alpha$, and $PLC$.
  • The equation for $C_a$ depends on variables $G_\alpha$, and $C_a$.

Thus one can detect a branching point at variable $G_\alpha$: there is one path leading to variable $PLC$ and one leading to variable $C_a$. The rate equations of $PLC$ and Ca do not directly depend on each other, only indirectly via $G_\alpha$.

Decomposition: Both branches of variable Gα are mathematically built in the same way. We can therefore decompose the model and only look at the subsystem $G_\alpha/PLC$. This is achieved by substituting a constant value for the concentration of calcium in the equation for $G_\alpha$ (e.g. $C_a = 0$). Surprisingly, the subsystem a-b still oscillates producing large amplitude spikes. These are reminiscent of the large spike component in the bursting.

image number 1

Equally if the concentration of $PLC$ in the equation for a is set equal to zero, the subsystem $a / C_a$ oscillates producing small amplitude fast oscillations. These resemble the secondary small-amplitude oscillations during the bursts.

image number 2

We conclude that the bursting oscillation is a complex oscillatory process that results from the interaction of two simple oscillatory processes that are in competition with each other. The analysis of the core model thus leads to a deeper understanding of the phenomenon of bursting in this signalling pathway.

Mathematics: Decomposing a model

The technique of understanding complex behaviour from its constituent components (sub-networks) can be illustrated in n example.

An equation of the form:

dx/dt = x - xy
dy/dt =     xy - y

has two variables $x$ and $y$ and constitutes interactions between the variables dur to the appearance of $y$ in the rate equation for $x$ and the appearance of $x$ in the rate equation for $y$. It can be decomposed by replacing individual variables with constant numbers:

If $y = 0$, the the first equation reduces to:

dx/dt = x

which means there is no impact of y and variable x will grow at a rate proportional to its current value.

If $y = 1$, the the first equation reduces to:

dx/dt = x - x = 0

meaning that variable $x$ does not change over time.

In this way we can understand the behaviour of the first equation for different ranges of the size of $y$ (assuming $y$ is positive as in population models):

  • If $y<1$, according to the fist equation variable $x$ will grow over time.
  • If $y>1$, according to the fist equation variable $x$ will shrink over time.

A similar argument applies to the rate equation for y if $x = 0$.

Decoupling sets of interconnected variables reduces the complexity of a network and sometimes makes counter-intuitive behaviour comprehensible.

This is a general method to study subsystems of a network and is particular useful to understand networks that comprise sub-networks in a modular structure.

7. Group 1: Pulse Perturbation
7.1. Simulations

Pulses of length 1 time unit can either leave the bursting behaviour mostly unaffected or cause a phase delay.

image number 1
image number 2
image number 3

The reason for the different impact can be seen in the state space plots. Red is the part of the trajectory during which the perturbation is applied. The perturbation at $t=12$ essentially does not have a visible impact of the trajectory in state space because $G_\alpha$, is near zero during that phase. The trajectory stays on its original course. The perturbation at $t=14$ causes a significant deviation of the trajectory after which it needs to settle to the original loop before resuming its original course.

MATLAB®: changing parameters during a simulation

A model function that is set up for simulations using a MATLAB® integration function like ode45, ode23 etc. requires the definition of a variable that represents time:

function dydt=model_Calcium(time, variables, parameters)

(where “parameters” are optional).

The independent variable “time” (or whatever name is assigned to it) can be accessed to define changes of parameters during the simulation.

The simulation period is defined in the script that calls the model function “model_Calcium” (e.g. $time = [0,50]$). Then, within the function you can use the variable “time” to put time-dependent changes of parameters in place:

k_2 = 2.9;

if time>12 && time<13
k_2 = k_2*1,3;

This will generate the stimulation protocol in the first of the pulse perturbation figures above.

image number 4
image number 5

So far, perturbation was defined as a temporary change in parameter value. There is a different way of simulating a perturbation: a change of concentration of one of the variables. This can be implemented by running a simulation to a certain point, resetting one or more of the variables, and running a second simulation with the reset values as starting condition.

image number 6

Here, the calcium concentration is increased at time $t=20$ and the simulation continued. The apparently small perturbation leads to the suppression of one of the bursts, i.e. it has a strong impact on the dynamics over the next 20 time units. From time $t=40$ onwards the bursting resumes as before.

MATLAB®: perturbation by a change of variables

Code to run state space perturbation (rather than parameter perturbation)

% Run first simulation with any starting condition:
[time1, vars1] = ode45(@model_Calcium, t_range1, var_init, options, k_1, k_2);

% Save last values as new starting conditions
var_init = vars1(end, :);

% Apply perturbation ca_stim to the third variable
var_init(3) = var_init(3) + ca_stim;

% Run second simulation
[time2, vars2] = ode45(@model_Calcium, t_range2, var_init, options, k_1, k_2);


  • Model name is “model_Calcium”, adjust to actual model function file.
  • Simulations require definition of $t_{range1}$, $t_{range2}$, $var_{init}$, $options$, $k_1$ and $k_2$.
  • Perturbation $ca_{stim}$ needs numerical value. In the case of decrease consider a factor rather than subtraction to avoid negative concentrations.
  • Before plotting, time and variable arrays need to be concatenated for display of the complete simulation.
7.2. Conclusion

Single pulse stimulation is a way to probe the dynamics of a model. It is useful to test predicted responses to perturbations that can be realised experimentally. In the case of rhythmic dynamics, simulations show how the timing, duration and strength of a perturbation may lead to qualitatively different outcomes.

Note: Following a suggestion by one participant, the model was checked for the presence of bistability (two different dynamics coexist at the same set of parameters) but none was found with the current settings.

8. Group 2: Periodic Perturbation

Periodic perturbation was implemented as the periodic modulation of one of the model parameters. For instance, we can replace $k_2$ with: $$k_2 \times (1+amp \times \mathrm{Sin}(fre \times time)$$

With periodic perturbation, two extreme cases are important to consider:

  • Frequency of the perturbation very slow compared to the processes that define the rate equations. In this case the observed dynamics is similar to what would be observed if the parameter was kept constant. Note however, how the permanent change in c) impacts on the number of secondary oscillations and the amplitude of the simple oscillations.
image number 1
  • Frequency of the perturbation very fast compared to the processes that define the rate equations. In this case the observed dynamics is similar to what would be observed if the parameter was kept constant at the mean of the perturbation. The model essentially suppresses the fast changes of the perturbation.

image number 2

At intermediate frequencies a complex interplay between the perturbation and the intrinsic dynamics takes place.

8.1. Impact of Amplitude

For constant frequency, resulting dynamics depends on the amplitude of the perturbation.

image number 3

a) There is regular bursting without perturbation. b) $amp=0.1$ has little impact on the bursting but impacts on the frequency such that a burst last for three perturbation periods. c) $amp=0.2$ disrupts the bursting structure and causes irregular bursting output. d) $amp=0.5$ causes spike-like responses which are interrupted by irregular “double spikes”. It appears that the resonse dynamics tries to follow the period of the perturbation.

The state space representation shows the impact of the perturbation on all three variables. In the illustrated case, the dynamics goes from regular to irregular, and there are three different return paths of calcium to base level, where the purple line hits the $PLC$ axis.

image number 4
8.2. Impact of Frequency

At constant amplitude, frequency has a huge impact on the both frequency and waveform of the response dynamics.

image number 5

Depending on the exact value, of the frequency, periodic perturbation can enhance the secondary oscillation of the bursting (top), induce new spiking patterns (centre), or substantially suppress the calcium signal (bottom). Unperturbed state as in Figure 12.

A scan of the frequency reveals that the response dynamics is extremely complex and no simple amplitude or frequency dependency on the perturbation frequency is to be expected according to this model.

image number 6

Here, multiple periods per stimulus frequency means that the response waveform is complex.

Given the range of time scales on which fluctuations in e.g. agonist concentrations occur we might expect typical calcium transients to be rather irregular in vivo. Conversely, downstream there might need to be particular mechanisms to deal such irregular unpredictable calcium transients when decoding this complex information.

9. Group 3: Signal decoding

The model of calcium dynamics can be extended to include a downstream process that depends on the calcium concentration and thereby acts as a decoder of the secondary messenger signalling. Following a study by Larsen et al the following equation was added to the model:

$$rate(EZ) = k_{act} \times \frac{C_{a}^P}{Km^P + C_{a}^P} - k_{inact} \times EZ$$

where EZ represents the activity of a calcium dependent enzyme.

$$EZ_{initial} = 0.01\\ k_{act} = 3\\ k_{inact} = 0.2\\ Km = 0.24\\ P = 1,4$$

This is a combination of one activation and one inactivation term. As expected, for a given presence of calcium, the output of this rate equation will be higher if $k_{act}$ increases and lower if $k_{inact}$ increases. What is not intuitive is to predict the response of this seemingly simple function to different temporal profiles of calcium.

9.1. Decoding of calcium profiles

Here is a graph of the average calcium concentration and the periods between peaks of the calcium signal. The range of parameter $k_2$ was chosen such that it contains non-oscillatory, simple spiking, complex bursting and overstimulating dynamics. Overall, calcium dynamics increases as $k_2$ increase, although the slope varies, being small initially and steepest for values above 3 for $k_2$.

The downstream response shows a remarkable fragmentation. Average enzyme activity remains near zero for low non-oscillatory calcium concentrations. Within the oscillatory region, enzyme activity increases near-linearly with $k_2$. Finally, there is a sudden switch to a very high activity level in the over-stimulation condition ($k_2>3$). The three distinct categories of response are enhanced for the cases where the Hill coefficient is larger than 1, i.e. in the case of cooperativity.

image number 1
image number 2
9.2. Dependence on $k_{act} / K_m$ ratio

The ratio $k_{act} / K_m$ determines the time scale of the activation of the calcium dependent enzyme activity. This time scale has profound consequences on the decoding process. Slow activation implies that repeated stimulation by periodic calcium variation can be summed up (integrated). The figure shows how for the case of small $k_{act}$ and $K_m$ value there is a steady increase to a high continuous level of enzymatic activity. In the case of large $k_{act}$ and $K_m$ values (but at same ratio) there is a transient rise followed by a decay to very small continuous levels. This is just one example of how two calcium dependent enzymes with different kinetic constants might pick up the same calcium signal dramatically different and thus lead to qualitatively different downstream responses.

image number 3
9.3. Dependence on $k_{inact}$

To understand this behaviour we can take a look at the function that depends on the values of these two parameters: $$f(k_{act},Km) = k_{act} \times Ca^P / (Km^P + Ca^P)$$. Here, the calcium concentration is assumed to be constant. The surface plot allows us to find combinations of parameters that give either a high or a low output for a given value of the calcium concentration.

% Figure plotted with:
ezsurf( 'X*1/(Y^4 + 1)', [0 10 0 3] )
image number 4
9.4. Dependence on $k_{inact}$

The dependence of the enzyme activity response on to changes in parameter $k_{inact}$ shows how the reaction rates define the filtering properties. For small values inactivation is slow and the average activity is high. Details of the fast oscillations in the calcium signal are smoothed out. At high values inactivation is fast and the average activity is low. However, the details of the fast calcium oscillation become apparent. Also the response now resembles a switch, enzyme activity is either turned “on” or “off”.

image number 5
9.5. Integration of fast oscillations

A curious effect can be seen when the nature of the calcium oscillation changes.

image number 6

The large-amplitude bursting leads to a period change in activity with an average of about 7 in the enzyme activity ($k_1=0.3$). The much smaller but faster calcium oscillations at $k_1=0.5$ counter-intuitively lead to accumulation of the activity to levels about twice as high.

10. Summary

The observation of two types of oscillatory responses in cytosolic calcium leads to a number of important insights into the dynamic complexities of living systems. This was studied using model equations based on a few key assumptions about molecular interactions. Basic mathematics and a number of programming techniques helped to illuminate the consequences of the assumptions in terms of expected biological features.

Focusing on a single example demonstrated that within a fairly short time one can get to surprisingly rich outputs from a computational model and formulate a large number of hypothesis for further experimental testing.

The SysMIC team would like to thank all participants for their engagement and contributions. Not all that was done could be included in this summary report but we have tried to capture the spirit of the work done at the Lane End Conference Centre.

Most of all we hope that the Workshop provided inspiration for further study of mathematics and computation in the service of biological research.

11. References
  • U. Kummer, L.F. Olsen, C.J. Dixon, A.K. Green, E. Bornberg-Bauer, G. Baier, Switching from Simple to Complex Oscillations in Calcium Signaling. Biophys J 79, 1188–1195 (2000).
  • A.Z. Larsen, L.F. Olsen, U. Kummer, On the encoding and decoding of calcium signals in hepatocytes. Biophys Chem 107, 83-99 (2004).
12. List of Parameter Values
$k_1$see text
$k_2$see text