Modeling calcium signaling by Cellular Automata simulation incorporating endocrine regulation and trafficking in various types of receptors

Calcium signaling plays important physiological roles which range widely from activation of muscle, synaptic transmission, cell movement, fertilization, growth in the cell population, learning process and retention of memory, to saliva secretion. Calcium also plays a crucial role in biochemical processes such as enzyme activity regulation, ion channels permeability, including activity of ion pumps, and different parts of the cytoskeleton. Signaling is initiated when the cell is stimulated, which results in the release of calcium ions (Ca2+) from intracellular stores. Many cell surface receptors, including G protein-coupled receptors, stimulate the formation of IP3 that diffuses to the endoplasmic reticulum, binds to its receptor, leading to the release of Ca2+ from the endoplasmic reticulum. Here, we construct a Cellular Automata model of signal transduction pathway mediated by calcium sensing receptors. The mechanism of receptor trafficking and dimerization is taken into account to exert positive feedback effects on the binding affinity of the cognate receptors. Difference equations are utilized to update the levels of calcium in the intracellular and extracellular compartments, incorporating endocrine regulation through the parathyroid hormone. Time series of calcium concentrations are obtained to investigate the effects of different fractions of healthy (functional) and defective receptors as well as dimerization process.


Introduction
Ever since the calcium-sensing receptor (CaSR) was identified in , a great deal of discoveries have been made regarding the role of calcium metabolism in health and disease []. Specifically, in being able to establish the role of calcium in the mechanism which regulates the secretion of parathyroid hormone (PTH), it has become possible to elucidate the genetic basis of numerous clinical disorders. According to Goodman [], the CaSR mediated signal transduction has now been pinpointed as an important signaling pathway in several types of tissues such as kidney, pancreas, gastrointestinal tract, placenta, and even brain tissues. It is now recognized that the extent to which the CaSR is activated modifies cell proliferation, cell differentiation, and apoptosis. In fact, CaSR affects the function of several key components of the parathyroid gland since disturbances in this signal-transduction pathway is believed to significantly contribute to the proliferation of parathyroid cells [].
In the parathyroid tissue, CaSRs are the detectors utilized by the parathyroid cells to sense variations in the extracellular calcium or blood ionized calcium concentration, consequently modulating parathyroid hormone secretion, so that serum calcium levels are maintained within a narrow physiologic range. When blood ionized calcium concentration is high, it binds to the receptors, which inhibits PTH secretion. When blood ionized calcium concentration is low, the receptor is unbound and inhibition on PTH secretion is lifted, resulting in its release.
In , Rattanakul and Lenbury [] proposed a stochastic Cellular Automata (CA) model of the signal transduction pathway modulated by the calcium sensing (CaS) receptors in order to study how different fractions of healthy (functional) and defective receptors modulate serum calcium concentration. In their model, receptor trafficking is assumed to take place on the cell surface so that dimerization and oligomerization occur, increasing receptor binding affinity. Comparing the time courses of bound receptors of various types with receptor trafficking with the case where there is no trafficking, it was found in [] that trafficking apparently facilitates receptor binding. The intracellular calcium fluctuates at a markedly higher level when receptor trafficking is allowed in the model. However, their model of calcium dynamics did not take into account the role of PTH in the feedback control mechanism and could not simulate oscillatory behavior in the extracellular calcium concentration, though this is often observed in experimental data [, ]. In this paper, we therefore propose an improved model system of difference equations to update the calcium concentrations at each time step that incorporates the effects related to PTH secretion. We also modified the probability expressions and parameter values used to update different types of receptors to reflect more accurately the stimulating effect and inhibitory effect of extracellular and cytosolic calcium, respectively. The CA model coupled with the model of calcium dynamics tracks, at time t, the numbers of in-activated receptors (N t ), free active receptors which are healthy (H t ), being the unbound receptors which function normally and are unbound, faulty active unbound receptors of type , F t (representing the free faulty receptors which possess lower binding affinity than H t and are not as long-lived), active defective unbound receptors F t (being very transient and short-lived, dissociating immediately after binding), ligand-bound healthy receptors (BH t ), bound defective receptors of type  (B t ), and bound defective receptors of type  (B t ).
Our model is shown to be capable of exhibiting oscillatory behavior in calcium concentrations. The time courses of different types of CaSR and calcium levels in the trafficking case are compared with those in the case where trafficking of receptors is not incorporated. As explained by Maurya and Subramaniam [], the profile of cytosolic calcium dynamics constitutes a valuable indicator of physiological events occurring at the cellular level. It also provides physicians with the means to measure the extent to which cells respond to external stimuli.

System model of calcium dynamics
In parathyroid cells, increases in extracellular calcium concentration activate the CaSR and diminish PTH release. In contrast, decreases in extracellular calcium concentration inactivate the CaSR and trigger the release of PTH that has been sequestered in secretory granules within the parathyroid cell. Activation of the CaSR induces transient increases in cytosolic calcium concentration owing to the mobilization of calcium ions from intracellular stores. Calcium ions are generally held in reserve within the endoplasmic reticulum (ER) until extracellular signaling causes the release of intracellular calcium into the cytosol [].
Rapid increase in cytosolic calcium transients is associated with decreased hormone PTH release. When extracellular fluid calcium is high, calcium binds to the receptor and this inhibits PTH secretion.
According to Shoback et al. [], 'the parathyroid cell is unusual because of the fact that low extracellular Ca + concentrations stimulate, while high Ca + concentrations inhibit, parathyroid hormone (PTH) release. ' Furthermore, Mayer and Hurst [] discovered from their experiment with anesthetized calves that, when plasma calcium concentrations became greater than . mg/ ml, a basal, non-diminishing PTH secretion rate of . ng/kg/min was maintained even if they induced a very high level of calcium. Small changes in PTH secretion rate were observed with changing plasma calcium if the calcium level stayed within the range between  and . mg/ ml. However, below  mg/ ml, a small decrease in plasma calcium concentration induced a marked increase in PTH secretion rate. According to Mayer and Hurst [], a maximum secretion rate of approximately . ng/kg/min is achieved at a plasma calcium concentration of about . mg/ ml.
Thus, we write the following difference equations to update cytosolic calcium concentration I t and extracellular fluid calcium concentration E t above their respective basal levels and PTH level X t at each time step t.
where H(·) denotes the Heaviside step function, ε is the rate at which intracellular calcium diffuses out of the compartment [], s is the scaling factor between extracellular ↔ intercellular compartments, δ  is the zero order intracellular calcium production rate, α is the intracellular calcium production rate constant arising from the heathy receptors, β and γ are those resulting from defective receptors of types  and , respectively. K F is the binding constant and K dF is the saturation constant in the fura binding process. K B is the buffering constant and K dB is the saturation constant in the calcium buffering process. These functions for calcium buffering and binding have been written based on the work of Blumenfeld et al. []. In (), ω is the rate at which calcium is removed from the extracellular compartment [], and the last term accounts for the rate at which extracellular calcium changes with the PTH level with rate constant λ. I b and E b are the basal values of the intracellular and extracellular calcium concentrations, respectively, and t stands for one time step.
so that the change in X t is zero once X t reaches its equilibrium value X e . When extracellular calcium drops below E M and where η  is the zero order removal rate of PTH by natural means.

Cellular Automata model of calcium signaling
According to Berto and Tagliabue [], CA models are discrete in space and time, consisting of a countable family of uniform, simple units, typically called either 'cells' or 'atoms' . At each time step, each cell is in a state within a finite set of states. At discrete time steps, the cells are updated according to specified transition rules and advance in parallel, depending on the status of neighborhood cells. Cellular Automata modeling is recognized as being flexible and offers researchers a better alternative when attempting to model a system involving multiple components.
neighborhood, H t may change into HB t and form a dimer if neighborhood, but there is one such receptor in its distant neighborhood, it can move closer to a randomized distant bound receptor, if there are more than one, and change A CA model was constructed for the calcium signaling process mediated by multiple types of receptors in [] and [], hence will not be described here in great detail. Here, we give, for completeness, the CA updating rules in Table , which follow those in [] and [], but with different probability expressions and different parameter values. This leads us to simulation results different from those presented in the previous work because the signal amplification and inhibitory effects are more properly incorporated as follows.
In Rule , the probability ρ  at which a non-activated receptor becomes activated is now which increases with increasing E t . In Rule , the probability ρ  , at which healthy receptors bind to calcium, increases with increasing extracellular calcium concentration, and also increases with increasing second messenger I t according to This is to take into account the amplification effect of the secondary messenger directly on the signaling strength.
In Rule , the probability ρ  at which the bound healthy receptors becoming unbound decreases as I t increases according to ρ  = e -k - e -k  I t and similarly for ρ  . This is to take into account the amplification effect of cytosolic calcium on binding affinity of the receptors. The negative feedback of cytosolic calcium is taken into account through its regulation by PTH in the updating equations ()-() that model calcium dynamics involving PTH concentration.
The simulation, performed in MATLAB, starts with the following initial conditions. E  denotes the initial number of empty units, N  is the initial number of non-active receptors, H  is the initial number of healthy active free receptors, whereas F  denotes the initial number of defective free receptors of type , F  is the initial number of defective free receptors of type , BH  is the initial number of bound healthy receptors, B  is the initial number of defective bound receptors of type , and B  denotes the initial number of defective bound receptors of type , the values of all of which are given in Table . According to Bajzer et al. [], at the steady state, the number of receptors at the surface of a cell is given by R  = V r k t , where V r is the zero order rate of insertion of receptors into the membrane and k t is the turnover rate constant of ligand-free receptors. Using the values of V r and k t found in [], we are led to R  ∼ = , receptors. Also, as explained  Table . All non-activated receptors are first updated at random. Next, all ligand-bound receptors are updated randomly, after which the rest of the cells are updated. At each time step, t =  s, we randomize a ρ,  ≤ ρ ≤ . In our Cellular Automata model, the cell membrane is represented by a square lattice having the size L × L, using Moore's neighborhood with a neighborhood of range .
For the updating rules above, we use parametric values estimated in [] and [] based on available literature. However, in the modified expressions, the parameters have been estimated to be in the feasible ranges reported in various sources so that the simulated time courses conform with those shown in the literature. These values are given in Table .
After all receptors have been updated at each time step, calcium and PTH concentrations are updated using the difference equations in ()-(). To explain more clearly how PTH feedback mechanism makes itself felt in our model, we show in Figure , the diagram of the interactive feedback loops at work in our simulations. Starting with the level of extracellular calcium, E c = E t + E b , equation () has been written as explained in the earlier section so that if E c increases in the range E m < E c < E M , the level of PTH will decrease. On the other hand, if E c decreases within this range, PTH level will increase. External to this range PTH remains more or less constant as has been observed experimentally.
The increase in PTH increases the rate of production of E t according to the last term in equation (). The variation in E t impacts on the chances of receptors being activated or being bound through the probability expressions ρ  and ρ  , respectively. The numbers   In (a), the receptor will move to the location (i -1, j + 1) if a bound receptor is found at one of the locations shown surrounded by double boxes on the upper left hand corner . In (b), the receptor will move to the location (i + 1, j + 1) if a bound receptor is found at one of the locations shown surrounded by double boxes in the upper right hand corner. In (c), the receptor will move to (i + 1, j -1) if a bound receptor (doubly boxed) is found in the lower right hand corner. In (d), it moves to (i -1, j -1) if a bound receptor (doubly boxed) is found in the lower left hand corner. Finally, in (e) the receptor will move either vertically or horizontally if a bound receptor is found vertically above (or below) it, or to its left (or right), respectively. of bound receptors of various types in turns lead to different levels of extracellular and intracellular calcium according to equations () and (), respectively. Thus, excessively high or low calcium level is brought back to a normal level by such PTH feedback mechanism.
We refer the readers to [] and [] for the detail of the manner in which the receptors are assumed to move toward neighboring receptors to form dimers or oligomers, after which their binding affinity is improved. Here, we show the diagram in Figure  to explain how a receptor moves closer to its neighboring receptors that are bound to form dimers or oligomers. If there are more than one bound receptor in its neighborhood, one of such receptors is chosen at random to be the one approached.

Simulation results and discussion
We discovered that the simulated time series from all simulation runs are not noticeably different from each other, and therefore, we show in the following figures the average over  runs for each construct. In Figure  the simulated time series of extracellular calcium concentration is shown, comparing the case when trafficking is not incorporated into our CA model (ρ  = .), corresponding to the dashed curve, to the case that receptor trafficking is incorporated, which is the solid curve. We see that both curves exhibit damped oscillations, while the receptor trafficking has the expected effect of significantly raising the extracellular fluid calcium concentration. Trafficking facilitates dimerization, which increases binding affinity of the receptors as has been reported by Wu et al. []. It is necessary for receptors such as the vascular endothelial growth factor to form dimers to become activated. Before the receptor monomers bind with the growth factor, they are   Tables 1 and 3. observed to move over the cell membrane, during which time they are called 'singletons' or take the form of pre-associated dimers which are ligand independent []. Figure  shows the simulated time series of intracellular calcium which also exhibits oscillations which dampen to a higher equilibrium level in the trafficking case. These graphs resemble many reported experimental data [, , -]. Figure  shows the simulated time courses of various types of bound receptors.
According to Linderman [], such a formation of clusters through dimerization influences the spatial organization of receptors across the cell membrane and, consequently, plays a crucial role in the signaling process. Evidence suggests that dimerization is used by cells to regulate the flow of information through the signal transduction pathways. Our model illustrates this modulation by dimerization formation and clustering, as seen in Figure , which shows a significant hike in the bound receptors, free or bound, when trafficking is taken into account.
In the study of type  Diabetes Mellitus (TDM), the disease is believed to arise from cognate receptor's resistance to insulin, a hormone secreted by the pancreatic cells. However, recent findings [] have suggested that the secretion of insulin by the beta cells is impaired in these TDM patients. Possibly, defective signaling leads to abnormal adjustments of calcium levels in circulation [].
In addition, Vauquelin [] investigated the slow dissociation of receptors and how that facilitates long-lasting protection of in vivo receptors. According to Vauquelin [], evidence has frequently been presented to support the proposal that the better antagonists are able to form complexes with their designated receptors which dissociate slowly, the more they are apt to give rise to longer receptor preservation. In this respect, we can see in our simulated graph of healthy bound receptors ( Figure (d)) here that its level is very much higher than the no-trafficking case. This is because the receptors in this figure are healthy and therefore slow in dissociation. Moreover, trafficking felicitates healthy receptors' binding with Ca + and, consequently, can contribute to long lasting actions on calcium kinetics (first term on the right of equation ()).
We also observe that, while the number of healthy bound receptors and type  defective receptors eventually plateau out to stable levels, the level of type  bound faulty receptors keeps on increasing, or it might take quite a bit longer to level off to a steady state level. This is because the bound receptors dissociate immediately the time step right after their binding to plasma ionized calcium since they are faulty of type . On unbinding, they can convert back into inactivated receptors or faulty receptors which are free to bind again and again very quickly, while the other types of receptors are more stable once bound to calcium. The pool of these types of free receptors is exhausted very quickly since there is insufficient time for it to be replenished. We can see that in the long run, the 'bad' receptors may rear up to create subsequent abnormalities in the functioning of several systems that rely on the proper functioning of calcium signaling, such as insulin secretion mechanism, proliferation of parathyroid cells and the development of parathyroid gland hyperplasia in chronic renal failure. The simulation runs shown here have been carried out with parameter values in the average ranges to validate that our CA model coupled with the calcium updating equations is able to exhibit dynamic behavior observed in measurements of normal individuals. Experimenting with other sets of parameters, we have been able to observe other time courses that might be found in an ailing patient. This allows us to differentiate health from sickness and identify the ranges of physical parameters which portend illness that could assist physicians in making better diagnosis and treatment decisions.

Conclusion
Our CA model combined with calcium updating difference equations is capable of simulating the dynamics of different types of receptors together with the concentration of intracellular calcium, believed to play a significant role in intercellular signaling [], and takes part in several important cellular functions that control many protein activities.
Previous modeling efforts have not incorporated the dynamics of PTH by closely observing the interactions between the extracellular ionized calcium with the fluctuations of PTH. Thanks to the detailed report by Mayer and Hurst [] on how the secreted PTH is affected by the changes in the calcium level, we have been able to incorporate more precisely the PTH compartment into the difference equation model of calcium levels, making use of the Heaviside step function.
Dynamic changes in cytosolic calcium and parathyroid hormone profiles provide us with a powerful tool to identify cellular activities since they constitute essential quantitative measures of the extent to which cells respond to external stimuli. Thus, our model and its simulations are meant to provide deeper understanding and form a basis for further studies to investigate how the function or dysfunction of this system impacts other intricate mechanisms in the human body, such as the connection of calcium signaling system to bone remodeling process and other related mechanisms, which is a subject for future study.