Computational Modelling of Nonlinear Calcium Waves
aa r X i v : . [ n li n . PS ] M a r Computational Modelling of Nonlinear Calcium Waves
Xin-She YangDepartment of Engineering, University of CambridgeTrumpington Street, Cambridge CB2 1PZ
Abstract
The calcium transport in biological systems is modelled as a reaction-diffusion process. Non-linear calcium waves are then simulated using a stochastic cellular automaton whose rules arederived from the corresponding coupled partial differential equations. Numerical simulationsshow self-organized criticality in the complex calcium waves and patterns. Both the stochasticcellular automaton approach and the equation-based simulations can predict the characteristicsof calcium waves and complex pattern formation. The implication of locality of calcium distri-bution with positional information in biological systems is also discussed.
Keywords: calcium transport, stochastic cellular automata, complex system, nonlinear waves,pattern formation, self-organized criticality.
PACS Numbers:
Citation detail:
X. S. Yang, Computational modelling of nonlinear calcium waves,
Applied Math-ematical Modelling , (2), 200-208 (2006). Many processes in living organisms such as muscle mechanics, cardiac electrophysiology, adaptationin photo-receptors, and cytoplasm functions involve the calcium ion transport and its physiologi-cal functions[1-6]. However, the exact function of Ca oscillations and transport is only partiallyunderstood although it is believed that they involve in the intracellular communications and syn-chronization in the response to a local stimulus [13-15]. Even in the simplest one-dimensional case,the proper modelling requires many simplifications and assumptions. The nonlinearity and cross-coupling in the transport mechanisms usually lead to intractable governing equations. Even withcertain approximations and simplicity, the mathematical models still lead to nonlinear reaction-diffusion equations if the essence of the calcium ion activities is included.Nonlinear reaction-diffusion systems can exhibit complex pattern formation [9-12]. The nonlinearsystem can be simulated using finite difference or finite element methods, or the alternative cellularautomata [16]. However, existing implicit numerical solution schemes are not always robust undergeneral boundary conditions. Most of the earlier work have mainly concerned the one-dimensionalcase with piecewise linear models [10]. Meanwhile, the cellular automata [17,18] have successfullymodelled reaction-diffusion systems [16] with relative stable pattern formation. This provides apossibility of dealing with the difficult nonlinear problem of calcium oscillation and waves usingfinite-state cellular automata with rules derived from related partial differential equations.On the other hand, self-organized criticality has been found in many systems in nature [7,8] sinceits discovery by Bak [7] and his colleagues. Since calcium oscillation and waves a very complicatedphenomenon with the characteristics of nonlinear reaction-diffusion systems, it can be expected thatregular patterns under certain conditions. One natural question related to this is: Do the self-organized criticality exist in the complexity of the calcium transport? However, this is no existingliterature addressing this question in the context of calcium transport. This is partly because the1esearch on calcium activities and their mathematical modelling for biological and physiologicalprocesses is still at a very earlier stage [10,14-15].This paper therefore first extends the existing mathematical models for nonlinear Ca waves.Then, a new stochastic finite-state cellular automaton is then formulated to simulate the generalnonlinear calcium transport process and pattern formation. Based on the reaction-diffusion equa-tions, the stochastic CA will be formulated. The pattern formation of calcium ions with realisticparameters will be studied. The self-organized criticality will be tested in the complex patterns ofcalcium concentration. The positional pattern formation and its implication will briefly be discussed. In many physiological processes within cells, Ca plays an essential role in controlling cellular be-havior and functioning in the sense that calcium ions act as a signalling agent for a wide range ofcellular activities. Calcium signaling is mediated through oscillation in intracellular Ca concentra-tion. Calcium ions can bind to a vast number of proteins and enzymes, and the binding can initiatea series of reactions that ends in the formation of a chemical called inositol 1,4,5-trisphosphate (IP ).The diffusion of Ca and IP through the cell cytosol can induce further release of calcium ionsfrom stores in the endoplasmic reticulum (ER) through IP -sensitive channels. These channels aresensitive to calcium itself, with fast activation for lower concentrations and comparatively slowerinhibition, thus leading to the calcium-induced calcium release (CICR). Complex wave characteris-tics such as plane waves and spiral waves have been observed in experimental studies using Xenopus oocytes as pointed out by McKenzie and Sneyd [19]. However, the detail mechanisms underlyingthese oscillation and waves are only partially understood.There are several mathematical models for Ca wave propagation in the literature, and theseinclude the important models such as DeYoung and Keizer [2], Goldbeter et al. [20], Sneyd et al. [15]and McKenzie and Sneyd [19]. However, different cell types can results in different mathematicalmodels, although the fundamentals are quite similar. The nonlinearity in the mathematical modelscan lead to complex pattern formations and spiral waves [21-24]. The model we will use in this paperis an extended version of a two-pool process. In the two-pool model, the calcium Ca concentrationin the cytoplasm and the concentration in the Ca -sensitive pool satisfy the two-pool model [1].Although it can reproduce the oscillations and waves observed in Xenopus oocytes and generatingspiral waves in the higher dimensional situations. The parameters for our simulations are based onthe models developed by Atri et al [1] and McKenzie and Sneyd [19]. The functioning of IP in theoscillation and waves has been reviewed in details by Sneyd et al. [25], and th calcium signaling hasbeen reviewed by Falcke [26].We use the variable u ( x, y, t ) =[Ca ] for the concentration ( µM ) of calcium, v ( x, y, t ) for thefraction of the IP receptors that are active. The variable p =[IP ] can be treated as a bifurcationparameter for the reason discussed later. Then, the nonlinear model equations for calcium transport[10,13-15] can be written as ∂u∂t = D u ∇ u + J f − J p + J l , (1) λ ∂v∂t = D v ∇ v + g ( u, v ) , (2)where J f models the flux of calcium through the IP receptor. J p models the amount the Ca being pumped out of the cytoplasm back into the endoplasmic reticulum or out through the plasmamembrane. J l models the calcium leaking into the cell. We have J f = kv [ δ + (1 − δ ) uk + u ] ε ( p ) , J p = γu d k du + u d , (3) J l = α, g ( u, v ) = k m k m + u m − v, (4)2nd ε ( p ) = p n k n + p n , (5)where k , k , k , k u are constants and α, β, δ are parameters. In addition, D v = 0 is used in mostexisting models since most models base on the assumption that Ca instantaneously activates theIP receptor.As the number of IP receptors remains approximately constant, thus we may assume that p isfixed and subsequently it can be considered as a bifurcation parameter. In fact, some studies using p ∈ [0 . , . ε ( p ) is the fraction of IP receptors that have bound IP and increases as p increases. Then, the complete nonlinear modelequations become ∂u∂t = D u ∇ u + kw [ δ + (1 − δ ) uk + u ][ p n k n + p n ] − γu d k du + u d + α, (6) λ dvdt = k m k m + u m − v, (7)where w is a time-dependent inactivation variable, and n = 3 , d = m = 2. Some typical values of therelated parameters are k = 3 µM s − (M=Mol/L is the molar concentration), k = k = 0 . µM, k =0 . µM, k v = 1 µM, k u = 0 . µM, δ = 0 . , α = 0 . µM s − , λ = 0 . s , D u = 20 µm s − , γ =2 µM s − [10].By using the notation s ( u, v ) = J f − J p + J l and g ( u, v ) = k m / ( k m + u m ) − v , we can obtain thesteady state solution ( u , v ). From equations (1) and (2) together with equations (6) and (7), it isstraightforward to check that they satisfy the Turing instability conditions Q = ( s u s v g u g v ) , (8)which requires tr( Q ) = s u + g v < , Det( Q ) = s u g v − f v g u > . (9)The model equations for intracellular calcium waves can be generally written in a system ofreaction-diffusion equations in the form φ t = D ∇ φ + f ( φ ) , φ = [ u v ] T , (10)where D = diag( D u , f ( φ ) is a general nonlinear function coupling the different compo-nents u and v . This nonlinear reaction-diffusion system can be solved using conventional numericalmethod or stochastic cellular automata. Conventionally, reaction-diffusion systems can be solved numerically using finite difference or finiteelement methods. Alternatively, the nonlinear systems can be simulated by using cellular automata[16]. The cellular automata for reaction-diffusion systems can be formulated to correspond to thesolutions of the related partial differential equations in a qualitative manner [4] or quantitativemanner [5]. The macroscopic approach via cellular automata as demonstrated by [16] is always moreefficient than explicit numerical method and can be more efficient than better numerical techniquein many cases.The discretization of equation (10) can be written as φ n +1 i,j − φ ni,j ∆ t = D [ φ ni +1 ,j − φ ni,j + φ ni − ,j (∆ x ) + φ ni,j +1 − φ ni,j + φ ni,j − (∆ y ) ] + f ( φ ni,j ) . (11)By taking ∆ t = ∆ x = ∆ y = 1, this becomes φ n +1 i,j = D [( φ ni +1 ,j + φ ni − ,j + φ ni,j +1 + φ ni,j − ) − φ ni,j ] + f ( φ ni,j ) , (12)3 Figure 1: Pattern formation of calcium concentrations (spiral waves on the left and rings on theright) using values given in section 2.1which can be written as a generic form φ n +1 i,j = r X k,l = − r c k,l φ ni + k,j + l + f ( φ ni,j ) , (13)where the summation is over the 2 r + 1 neighborhood. The coefficients c k are determined from thediscretization of the governing equations, and for this special case, c − , = c +1 , = c , − = c , +1 = D, c , = − D . In fact, this is equivalent to a finite-state cellular automaton with a transition rule G = [ g ij ] ( i, j = 1 , , ..., N ) from one state Φ n = [ φ nij ] ( i, j = 1 , , ..., N ) at time level n to a newstate Φ n +1 = [ φ n +1 ij ] ( i, j = 1 , , ..., N ) at time level n + 1. That is, G : Φ n Φ n +1 , g ij : φ nij φ n +1 ij . (14)The basic assumption for the rule inference for stochastic automata is that the function g ij ( φ nij ) = φ n +1 ij , V ≤ Γ( φ nij , f, r, N ) , (15)where Γ is a function with a range of [0 , V is a random variable. At every time step, arandom number V is generated for each automaton ( i, j ). The new state will only be updated if thegenerated random number is greater than Γ, otherwise, it remains unchanged. Following a similarderivation for ordinary differential equation [9], the probability function Γ can be written asΓ( φ nij , f, r, N ) = Γ( e − f ) = a + be − f , (16)where a and b are coefficients depends on the number of finite states and other factors such asthe accuracy of the approximation to the partial differential equations. The parameters of thecontinuum reaction-diffusion model shall be calibrated to fit the results given by the stochasticcellular automaton model using a least-squared procedure so as to get the related accurate transitionrules. By using the approach of stochastic cellular automata with finite number of states, we can now sim-ulate the reaction-diffusion calcium transport and nonlinear calcium waves. Numerical simulationsare carried out on an N × N lattice in 2-D setting, and usually, N ≥
40, or up to 1500. The number4f states is taken to be 256 in the present simulations. Different simulations with different latticesize are compared to ensure the simulated results are independent of the lattice size and time steps.In the rest of the paper, we present some results related to calcium pattern formation, nonlinearcalcium waves and self-organized criticality of the complexity of reaction-diffusion systems.
From the initial random configuration, nonlinear reaction-diffusions can lead to complex patterns.Figure 1 shows two examples of the calcium concentrations evolving to the interesting patternsin 2-D configuration with 200 ×
200 cells. The parameters used in the calculations are k =3 µM s − (M=Mol/L is the molar concentration), k = k = 0 . µM, k = 0 . µM, k v = 1 µM, k u =0 . µM, δ = 0 . , α = 0 . µM s − , λ = 0 . s . Spiral waves are formed for values D u = 5 µm s − (Fig1a) while rings and ribbons are formed for higher value D u = 25 µm s − (Fig 1b). These patternsgradually evolve with time; however, the general characteristics of patterns only change slowly withtime.The simulations imply that patterns and structures formed by local transition rules are relativelystable. Our present results are consistent with the other well-known studies in pattern formationsuch as [11,12] for plants and sea shells by using nonlinear reaction-diffusion equations. In addi-tion, our present simulations suggest that patterns arise naturally from the local interactions eitherthrough rule-based/agent-based evolution in terms of relationships among the neighbourhood indi-viduals in stochastic cellular automata or partial differential equations in terms of system variablessuch as calcium concentrations. The initial configuration is not important. The rules of interac-tions or relationships between entities/individuals in the nearest neighbour are crucial factors thatresponsible for the behaviour of the system and pattern formations.In addition, the spatial pattern formation of calcium concentration can provide some positionalinformation of calcium distribution resulting from the reaction-diffusion transport among cells. Thismay have some important implication in the mechanism of calcium functionality in biological sys-tems. Although many factors such as the function of proteins, genetic information and enzymesaffect the intracellular transport of calcium, however, the calcium reaction-diffusion process itselfcould greatly contribute to the spatial distribution of calcium and thus be responsible to some extentfor its positional information and signalling in biological systems. If it is the case, the modelling ofcalcium transport can be beneficial to the understanding of formation of the spatial structure andpositional signalling coupled with the genetic and functional information in biological systems andphysiological mechanisms.The reaction-diffusion systems of calcium transport are complex. Nonlinear waves can ariseunder certain initial and boundary conditions. For a single source, the calcium concentration startsto expand and form nonlinear waves as shown in the figure on the left of Figure 2 which is asnapshot at t = 500. Numerical simulations imply that wave speed decreases with time as expect forany passive diffusion systems. In addition, the for a random configuration with periodic boundaryconditions, complex nonlinear wave pattern have observed in numerical simulations as demonstratedin the figure on the right (Fig. 2b) where calcium concentration varies with space and time. For a lattice N = 100 ×
100 with 256 states, the complexity of the cellular automata can be measuredby its entropy S is defined as S = − P i p i log p i , where p i is the probability of states i [1]. For afinite state automaton, p i can be approximated by the calcium concentration so that S = − X i u ij log u ij . (17)The variation of complexity or entropy of the stochastic cellular automata with time is shown inFigure 3a. It is clearly seen that the complexity varies significantly at the early stage of the patternformation process, then it gradually relaxes to the equilibrium at long time, indicating that thereaction-diffusion system is in a quasi-steady state.5igure 2: Nonlinear calcium waves for a single source (left) and multiple sources (right). e n t r o p y −2 −1 Calcium concentration N u m b er o f C e ll s Figure 3: Complexity and entropy variations at different time steps (left). Self-organized criticalityin calcium pattern formation leads to the exponent γ = 1 . ± .
02 (right).The complex pattern of calcium concentration can be measured by grouping or counting thenumber of cells for a given value of concentration. The results are plotted in Figure 3b. It is clearlyseen that there exists a power law in the distribution of the number of cells ( n ) versus discrete calciumconcentration ( u ). A least-square fitting of n ∝ u − γ , leads to the exponent of γ = 1 . ± .
02. Thisimplies the self-organized criticality in the complex calcium patterns.This result may have important implications to the calcium transport mechanism. Althoughthe detail of intracellular calcium oscillation and communication is not clear, it is likely that theintracellular and intracellular interactions mainly occur locally. Thus reaction-diffusion dominatesthe process without much contribution from the convection mechanism. This is consistent with thephysiological aspects of calcium transport and functioning [4,10].
The finite state stochastic cellular automata have been formulated to simulate the reaction-diffusionsystems of nonlinear calcium oscillations and waves using the transition rules derived from the partial6ifferential equations. By using the proper stochastic and transition rules between different states,the finite state automaton can simulate the complexity of calcium transport. Numerical experimentsshow that complex patterns can arise from the initial random configuration due to the local transitionrules between entities of certain nearest neighborhood and the details of initial configuration is notimportant. The power-law relationship between number of cells and calcium concentrations impliesself-organized criticality in the complex patterns. [1] Atri A., Amundson J., Clapham D., and Sneyd J., A single-pool model for intracellular calciumoscillations and waves in the Xenopus laevis oocyte,
Biophysical Journal , (1993) 1727-1739.[2] De Young G. W. and Keizer J., A single pool IP -receptor based model for agonist stimulatedCa oscillations, Proc. Natl. Acad. Sci. USA , (1992) 9895-9899.[3] Dupont G., Goldbeter A., Properties of intracellular Ca waves generated by a model basedon Ca -induced Ca release, Biophysical Journal , (1994) 2191-2204.[4] Ermentrout G. B., Edelstein-Keshet L., Cellular automata approaches to biological modelling, J Teor. Biol. , (1993) 97-133.[5] Gerhardt M., Schuster H., Tyson J.J., A cellular automaton model of excitable media includingcurvature and dispersion, Science , (1990) 1563-1566.[6] Goldbeter A., Biochemical Oscillations and Cellular Rhythms: the Molecular Bases of Periodicand Chaotic Behaviour , Cambridge University Press, Cambridge, (1996).[7] Bak P, Tang C and Wiesenfeld K, Self-organized criticality: an explanation of 1/f noise,
Phys.Rev. Lett. , (1987) 381-384.[8] Bak P, How nature works: the science of self-organized criticality , Springer-Verlag, New York,(1996).[9] Guinot V.,Modelling using stochastic, finite state cellular automata: rule inference from contin-uum model,
Appl. Math. Model. , (2002) 701-714.[10] Keener J., Sneyd J., Mathematical Physiology , Springer, (1998).[11] Murray, J. D.,
Mathematical Biology , Springer, New York, (1989).[12] Meinhardt H.,
Models of biological pattern formation , Academic Press, London, (1982).[13] Sanderson M. J., Charles A. C., Boitano S., and Dirksen E. R., Mechanisms and function ofintercellular calcium signaling,
Molecular and Cellular Endocrinology , (1994) 173-187.[14] Sneyd J., Charles A. C., and Sanderson M. J., A model for the propagation of intercellularcalcium waves, Am J Physiol. , (1994) 293-302.[15] Sneyd J., B Wetton, Charles A. C., and Sandeson M. J., Intercellular calcium waves mediatedby diffusion of ionositol trisphoshpate: a two-dimensional model, Am. J. Physiology CellPhysiology , (1995) 1537-1545.[16] Weimar J. R., Cellular automata for reaction-diffusion systems, Parallel computing , (1997)1699-1715.[17] Wolfram, S., Cellular automata as models of complexity, Nature , (1984) 419-424.[18] Wolfram, S., Cellular automata and complexity , Reading, Mass: Addison-Wesley (1994).719] McKenzie A. and Sneyd J., On the formation and breakup of spiral waves of calcium,
Inter-national Journal of Bifurcation and Chaos , (1998) 2003-2012.[20] Goldbeter A., Depont G., and Berridge M., Minimal model for signal-induced calcium oscil-lation and for their frequency encoding through protein phosphorylation, Proc. Natl. Acad.Sci. USA , (1990) 1461-1465.[21] Hunding A. and Ipsen M., Simulations of waves in calcium models with 3D spherical geometry, Mathematical Biosciences , (2003) 45-66.[22] Romeo M. M. and Jones C. K. R. T., The stability of travelling calcium pulses in a pnacreaticacinar cell, Physica D , (2003) 242-258.[23] Wilkins M. and Sneyd J., Intercellular spiral waves of calcium, J. Theor. Biol. , (1998)299-308.[24] Barkley D., Linear stability analysis of rotating spiral waves in excitable media, Phys. Rev.Lett. , (1992) 2090-2093.[25] Sneyd J., Falcke M, Dufour F. F., Fox C., A comparison of three models of the inositoltrisphosphate recptor, Prog. Biophys. Mol. Biol. , (2004) 121-140.[26] Falcke M, Reading the patterns in living cells – the physics of Ca signaling, Adv. Phys. ,53