Buffering capacity explains signal variation in symbiotic calcium oscillations.

Legumes form symbioses with rhizobial bacteria and arbuscular mycorrhizal fungi that aid plant nutrition. A critical component in the establishment of these symbioses is nuclear-localized calcium (Ca2+) oscillations. Different components on the nuclear envelope have been identified as being required for the generation of the Ca2+ oscillations. Among these an ion channel, Doesn't Make Infections1, is preferentially localized on the inner nuclear envelope and a Ca2+ ATPase is localized on both the inner and outer nuclear envelopes. Doesn't Make Infections1 is conserved across plants and has a weak but broad similarity to bacterial potassium channels. A possible role for this cation channel could be hyperpolarization of the nuclear envelope to counterbalance the charge caused by the influx of Ca2+ into the nucleus. Ca2+ channels and Ca2+ pumps are needed for the release and reuptake of Ca2+ from the internal store, which is hypothesized to be the nuclear envelope lumen and endoplasmic reticulum, but the release mechanism of Ca2+ remains to be identified and characterized. Here, we develop a mathematical model based on these components to describe the observed symbiotic Ca2+ oscillations. This model can recapitulate Ca2+ oscillations, and with the inclusion of Ca2+-binding proteins it offers a simple explanation for several previously unexplained phenomena. These include long periods of frequency variation, changes in spike shape, and the initiation and termination of oscillations. The model also predicts that an increase in buffering capacity in the nucleoplasm would cause a period of rapid oscillations. This phenomenon was observed experimentally by adding more of the inducing signal.

Plant growth is frequently limited by the essential nutrients nitrogen and phosphorus. Several plant species have established symbiotic relationships with microorganisms to overcome such limitations. In addition to the symbiotic relationship with arbuscular mycorrhizal fungi that many plants establish in order to secure their uptake of water, phosphorus, and other nutrients (Harrison, 2005;Parniske, 2008), legumes have developed interactions with bacteria called rhizobia, resulting in the fixation of atmospheric nitrogen within the plant root (Lhuissier et al., 2001;Gage, 2004;Riely et al., 2006).
Root symbioses initiate with signal exchanges in the soil. Plant signals are perceived by the symbionts, triggering the successive release of microbial signals. For rhizobia, the signal molecules are lipochitooligosaccharides termed Nod factors (Dénarié et al., 1996), and the release of lipochitooligosaccharides has also been found in the fungal interaction (Maillet et al., 2011). Upon receiving diffusible signals from the microsymbionts, the plant roots initiate developmental programs that lead to infection by rhizobia or arbuscular mycorrhizal fungi. Both programs employ the same signaling pathway with several components being common to both mycorrhizal and rhizobial interactions (Kistner and Parniske, 2002;Lima et al., 2009). In particular, both the symbioses show characteristic perinuclear and nucleoplasmic localized calcium (Ca 2+ ) oscillations, so-called Ca 2+ spiking Sieberer et al., 2009). It has been suggested that Ca 2+ is released from an internal store, most likely the nuclear lumen and associated endoplasmic reticulum (ER; Matzke et al., 2009), with targeted release in the nuclear region (Capoen et al., 2011).
Genetic screens in the model legume Medicago truncatula have identified several genes that are required for the plant in the establishment of both symbioses. Two of these, Doesn't Make Infections1 (DMI1) and DMI2, are essential for the induction of the Ca 2+ oscillations, yet the precise roles and mechanisms of these components remain to be determined. DMI2 codes for a plasma membrane receptor-like kinase (Endre et al., 2002;Stracke et al., 2002) that is required to facilitate further signal transduction in the cell (Bersoult et al., 2005). DMI1 is a cation channel located preferentially on the inner nuclear envelope Edwards et al., 2007;Riely et al., 2007;Charpentier et al., 2008;Capoen et al., 2011;Venkateshwaran et al., 2012). DMI3 encodes a calcium calmodulin-dependent protein kinase that interacts with downstream components and is thought to be the decoder of the Ca 2+ oscillations Mitra et al., 2004;Hayashi et al., 2010). Additional genes coding for three nucleoporins called NUP85, NUP133, and NENA are also required for Ca 2+ oscillations (Kanamori et al., 2006;Saito et al., 2007;Groth et al., 2010). The nuclear pore might be involved in trafficking secondary signals and/or ion channels to the inner nuclear membrane. These shared signaling components are collectively referred to as the common Sym pathway.
DMI1 plays a key role in the production of Ca 2+ oscillations, but its exact mechanism is still unknown. Orthologs of DMI1 have been found; in Lotus japonicus, they are called CASTOR and POLLUX (Charpentier et al., 2008), and in pea (Pisum sativum), SYM8 . CASTOR and POLLUX, as well as calcium calmodulin-dependent protein kinase, are highly conserved both in legumes and nonlegumes (Banba et al., 2008;Chen et al., 2009). This highlights the essential role of the Ca 2+ oscillations in mycorrhizal signaling.
DMI1 is not the channel responsible for the release of Ca 2+ (Charpentier et al., 2008;Parniske, 2008;Venkateshwaran et al., 2012) but probably influences the activity of Ca 2+ channels. This is similar to how some K + channels act in other plants and yeast (Peiter et al., 2007). Indeed, DMI1 is possibly a K + -permeable channel, based on the observation that POLLUX complements K + transport in yeast (Charpentier et al., 2008). In symbiosis, the mode of action of DMI1 could be to allow cations into the nuclear envelope and in that way counterbalance the transmembrane charge that would occur following the release of Ca 2+ into the nucleoplasm and cytoplasm. The cation channel could thus change the electrical potential across the nuclear membranes, affecting the opening of the voltageactivated Ca 2+ channels . This hypothesis is supported by a study reporting a membrane potential over the nuclear envelope in plants (Matzke and Matzke, 1986).
Pharmacological evidence and the characteristics of the Ca 2+ oscillations supports the involvement of Ca 2+ pumps and Ca 2+ channels (Engstrom et al., 2002). The pumps are needed to resequester Ca 2+ after each release event, actively transporting Ca 2+ against the concentration gradient using ATP. A recent study found a SERCA-type Ca 2+ ATPase, MCA8, that is located on the inner and outer nuclear envelope of M. truncatula and is required for the symbiotic Ca 2+ oscillations (Capoen et al., 2011). Such SERCA pumps are widely distributed on plant membranes, and the variation in their structure points to them being differentially regulated (Sze et al., 2000).
Ca 2+ channels release Ca 2+ from the store, and the mechanism of activating these Ca 2+ channels has been hypothesized to be voltage gated Oldroyd and Downie, 2008), but this remains to be verified experimentally. After release of Ca 2+ into the cytosol and nucleoplasm, buffers quickly bind to and remove these free ions due to their toxicity to the cell (Sanders et al., 2002). Buffers, i.e. molecules that can bind Ca 2+ , may play an important role in determining the nonlinear behavior of the oscillatory system for Ca 2+ signaling (Falcke, 2004). As numerous Ca 2+ buffers are present in cells, it is important to take their contribution into account. Such buffers can also include experimentally introduced dyes and Ca 2+ chelators.
In Capoen et al. (2011), we investigated the establishment and transmission of spatial waves across the nuclear envelope and demonstrated that the key components for Ca 2+ spiking reside on the inner and outer surface of the nuclear membrane. The computational framework we employed for this analysis makes a number of approximations in order to provide the computational efficiency required to perform spatiotemporal simulations. Here, a main focus is to understand the effect of buffers on the Ca 2+ oscillations.
In this article, we propose a mathematical model based on three key proteins; a Ca 2+ ATPase, a voltagegated Ca 2+ channel, and the cation channel DMI1. The model reproduces the symbiotic Ca 2+ oscillations, and we further demonstrate that Ca 2+ -binding proteins can explain initiation, termination, and experimentally observed variation in oscillation patterns. Furthermore, the model predicts that increases in buffering capacity can cause a period of rapid oscillations, and these were observed experimentally.

Assumptions Made to Model Ca 2+ Behavior
The model consists of three key proteins: a K + channel, a voltage-gated Ca 2+ channel, and a Ca 2+ ATPase. This is illustrated in Figure 1A, as well as these components' location on the inner nuclear membrane. A more detailed explanation of the nuclear localization is shown in Supplemental Figure S1. The K + channel represents the protein DMI1 (Venkateshwaran et al., 2012), which is modeled as a driver of changes in membrane voltage (Table I; Supplemental Fig. S5). The Ca 2+ channel is assumed to be voltage dependent based on studies supporting such channels in plants (Grygorczyk and Grygorczyk, 1998). The released Ca 2+ can feed back and bind to the K + channel, based on a study predicting partial gating of DMI1 by Ca 2+ . In animal systems, some Ca 2+ -gated K + channels are voltage gated, but the lack of voltage sensor motifs in the transmembrane domain of DMI1 suggests that changes in voltage would have no effect on the conductance of DMI1. The Ca 2+ pump is modeled to represent MCA8, a SERCA-type Ca 2+ ATPase (Capoen et al., 2011), and it resequesters the released Ca 2+ into the Ca 2+ store against the concentration gradient. In this model, DMI1 is presumed to be able to change the membrane voltage to provide a counterion current to the Ca 2+ movement.
The model assumes an inner nuclear envelope localization in M. truncatula root hair cells, based on the preferential localization of DMI1 (Capoen et al., 2011). The nuclear envelope lumen is contiguous with the ER and is assumed to be the Ca 2+ store, generating a concentration gradient across the inner nuclear envelope with high Ca 2+ levels inside (nuclear envelope lumen) and low Ca 2+ levels outside (nucleoplasm). For simplicity, we assumed an infinite Ca 2+ store, and there is only one pool of Ca 2+ in the model. That pool is the Ca 2+ in the nuclear envelope lumen and ER, which can be released into the nucleoplasm.

A Mathematical Model to Recapitulate Ca 2+ Oscillations
Our aim is to understand how the nuclear Ca 2+ concentration varies over time. We therefore chose to model the symbiotic Ca 2+ oscillations with a set of ordinary differential equations (ODEs). ODEs are a convenient way of describing time-dependent behavior. In addition, this framework allows us to build on published results of channel characteristics as a function of current and voltage.
The simple model consists of two ODEs (Table I). The first describes the changes in voltage over the inner nuclear membrane. The second describes the changes in free Ca 2+ inside the nucleus. A part of these ODEs is two currents: the Ca 2+ current through the Ca 2+ channel and the K + current through the K + channel. This model is electrical in nature and is similar to a model for action potentials in a pancreatic beta Figure 1. A, Schematic overview of the basic model. The three key components in the model are the K + channel, a voltagegated Ca 2+ channel, and a Ca 2+ pump depicted on the inner nuclear membrane. B, Experimental Ca 2+ oscillations using microinjection of the dyes OG and TR into a M. truncatula root hair cell. The trace has been detrended with the moving average method. C, Simulated Ca 2+ oscillations from the basic model shown in A, excluding buffers. D, Phase space diagram of Ca 2+ concentration and voltage of the oscillations shown in C. The system trajectory from an arbitrary initial condition is shown as a solid line making one full oscillation. The nullclines for the calcium concentration, c : = 0, and for the voltage, cell (Atwater et al., 1983;Chay and Keizer, 1983;Bertram et al., 2000;Fridlyand et al., 2003Fridlyand et al., , 2009Fridlyand et al., , 2010Bertram and Sherman, 2004;Tsaneva-Atanasova et al., 2006). Additional ODEs can be added to this model to account for changes in buffers that bind to Ca 2+ , and this is described by the term R i , where i represents an individual buffer species. In the simple model outlined above, R i is zero. When buffers are present in the model, we assumed they bind Ca 2+ in a 1:1 ratio. These buffers represent a perturbation to the basic system presented above. Perturbations such as changes in buffer capacity of the cell can influence the oscillations. Both of the model behaviors, oscillating or resting, are stable in the model; therefore, small perturbations only change the behavior quantitatively (Katok and Hasselblatt, 1995). For example, the introduction of buffers can lead to a change in oscillation frequency. However, when the change in buffer characteristics increases and they become a larger perturbation, then they can also change the time of transition between the two behaviors: oscillations and resting state.
The voltage-dependent Ca 2+ channel is modeled with the classical Hodgkin-Huxley-type gates (Izhikevich, 2007). The Ca 2+ channels conductance is voltage dependent, while the conductance of the K + channel is Ca 2+ dependent but not voltage dependent. The Ca 2+ pump is modeled as electrically neutral, with a constant pump rate, comparable to what has been described for many Ca 2+ ATPases (Thomas, 2009); therefore, in the model the Ca 2+ ATPase does not contribute to the changes in membrane potential. Four of the model parameters have been taken from previous measurements: the volume of the nucleus (Brière et al., 2006), the capacitance of the nuclear envelope (Grygorczyk and Grygorczyk, 1998), the resting potential of Ca 2+ (Petersen et al., 1998), and a scaling factor relating total Ca 2+ changes to changes in free Ca 2+ (Brière et al., 2006). The remaining parameters have been estimated, and details on parameter estimation methods can be found in "Materials and  Tables II and III and in Supplemental  Table S1 for the supplemental figures.
In dynamical systems theory, oscillations are called limit cycles and the resting states are called fixed points. We will use this terminology where appropriate. Table I. Equations for the oscillatory calcium model v is the voltage difference across the interior nuclear membrane, c is the concentration of Ca 2+ in the nucleus, and p i denotes the concentrations of N various buffers for i = 1,.,N. The association and dissociation constants for the i th buffer are labeled k þ i and k 2 i , and p tot i is the total concentration of the i th buffer. We assume binding to Ca 2+ occurs in a 1:1 ratio. The Ca 2+ current through the voltage-gated channel is described by I ca , and the K + current through the K + channel is described by I k . Definitions of all parameters are given in Tables II  and III

Implications of the Model
We first evaluated the initial model, lacking buffers, to test the behavior with a minimal set of components. This can provide some valuable insights and clues of how to extend the system to account for further observations. Even this simple setup produced oscillations that correspond well to those experimentally observed (Fig. 1, B and C), suggesting that these three components can create a self-sustaining oscillation. This was achieved by fitting the parameters such that the resulting Ca 2+ concentration profile from the model matched a single Nod factor-induced spike. The result is an oscillation, a stable limit cycle, with a period of approximately 85 s that oscillates between 130 and 630 nM Ca 2+ .
This oscillation is pictured in phase space in Figure  1D. The system also has three points where it does not oscillate, so-called fixed points, in this parameter range at (v,c) = {(55, 0), (54.975, 0.03185), (22.64, 0.158)}, each of which are unstable except for the point at (55, 0). If the concentration of Ca 2+ at the start of the simulation is below a threshold of approximately 22 nM Ca 2+ , then the system decays to this stable fixed point. The fixed points appear as circles in Figure 1D, while the nullclines for v and c appear in blue and green, respectively. The nullclines are the values for which the voltage or Ca 2+ do not change over time. Where nullclines intersect, neither of the system's variables change, in our case voltage and Ca 2+ ; thus, the system is at rest. In Figure 1D, notice the characteristic inverted "N" shape of the voltage nullcline and its single intersection with the Ca 2+ nullcline between the maximum and minimum of the "N." These features are often associated with oscillations in other models, such as the Morris-Lecar equations (Fall, 2002).
This simple model fails to capture several experimental observations. First, it cannot reproduce differences in spike shapes that are observed using different experimental techniques nor the variation in period seen between and within individual time series. Second, the model fails to mimic a common observation of faster frequency spiking at the start of oscillations. Third, the model will oscillate continuously and thus lacks a termination mechanism. Finally, the experimental data are not real Ca 2+ concentrations, but a relative measure of how much Ca 2+ is bound to the Ca 2+ reporter, an introduced buffer, so we should ideally fit the experimental data to Ca 2+ bound to buffers in the simulations. We therefore added Ca 2+ -binding buffers to the model and investigated their effect on oscillations.
Buffering Capacity: A Key Player in Shaping the Ca 2+ Signature?
We investigated the effects of changing the buffer characteristics on the dynamics of the system. Given that DMI1 binds Ca 2+ in our model, leading to changes in membrane voltage and the subsequent release of a more Ca 2+ , it is reasonable to consider that buffers binding to the released Ca 2+ might affect the oscillations.
When Ca 2+ levels oscillate in the nucleoplasm, the proportion of bound and unbound buffers will also oscillate, since more buffers will bind Ca 2+ when more Ca 2+ ions are available. The bound buffer could represent a Ca 2+ reporter, whose fluorescence indicates the level of Ca 2+ . However, we observed scenarios in our model where the oscillations of bound buffer do not precisely follow the Ca 2+ oscillations. At low offrates, k 2 , of Ca 2+ from the buffer, the oscillation of the bound buffer lags behind the Ca 2+ oscillations by several seconds (Supplemental Fig. S2A). However, as off-rates and dissociation constants, K d , are increased, the buffer pattern becomes more and more similar to the Ca 2+ spikes. Another scenario in which the buffer does not follow the Ca 2+ oscillations is when the basal level of Ca 2+ and the buffer concentration are higher than the amplitude of the Ca 2+ oscillations, and k 2 and K d of the buffer are small. In this case (Supplemental Fig. S2B), the concentration of bound buffer initially builds up to the level around which the bound buffer will oscillate. During this phase, a few high frequency Ca 2+ spikes are seen before the Ca 2+ oscillations settle into a periodic pattern with narrow peaks and long resting periods (about 200 s). However, for the current model with the parameters listed in Tables II and III, we found that the concentration of bound buffers accurately followed the Ca 2+ dynamics if there was more than one buffer in the system, as is demonstrated in Supplemental Figure S2C, so the model was extended to include two buffer species.
To measure Ca 2+ , we use two main experimental techniques of Ca 2+ indicators. The first is to inject fluorescent Ca 2+ -sensitive dyes into the cell, which increase their fluorescence when binding Ca 2+ (Rudolf et al., 2003). The second method is to transform the plants with Ca 2+ reporter genes that encode proteins that allow increases in Ca 2+ to be measured using fluorescent resonant energy transfer (Miyawaki et al., 1999). These two methods show differences in spike shapes as a result of different binding characteristics of the Ca 2+ indicators (Fig. 2, A and C). To allow the model to reproduce these experimental data, we optimized the parameters by matching model-generated spikes with different experimentally observed spikes (Fig. 2, B and D). The dissociation constant of one of the buffers was chosen to reflect the Ca 2+ indicators typically used in experiments, which is either Ca 2+ dyes or a Ca 2+ reporter. When the Ca 2+ dyes are used, both a sensitive dye and an insensitive dye are injected, and their ratio is taken. This decreases the impact of cytoplasmic streaming on the results. The Ca 2+ -sensitive dye Oregon Green (OG) has a K d of 0.17 mM, and the insensitive dye Texas Red (TR) a K d of 0.76 mM (Molecular Probes), while the Ca 2+ reporter Cameleon YC2.1 has a K d of 5.4 mM (Miyawaki et al., 1999). By changing the buffer characteristics to these known values in Figure 2, we were able to match both OG/TR and Cameleon YC2.1 spike shapes. Experimental traces show a significant variation in period both between and within time series, and by including buffers, this behavior could be captured. In general, large amounts of buffers with large dissociation and association rates cause an increase in the period (Fig. 3). The bifurcation diagrams in Supplemental Figures S3 and S4 summarize the impact of buffers on the oscillations.

Prediction and Experimental Validation of Rapid Ca 2+ Oscillations
Fast oscillatory transients are often observed in the initial phases of experimental Ca 2+ time series (Fig.  4A). To induce oscillations in M. truncatula root hairs, the Nod factor signal from M. truncatula's symbiont Sinorhizobium meliloti was used, henceforth referred to as Nod factor. After Nod factor application at a concentration of 1 nM, we observed rapid spiking in seven out of 13 cases. We simulate this phenomenon in the case of two buffers in Figure 4B. These periods of rapid spiking in the model simulation are caused by the presence of high levels of unbound buffers, and the rapid spiking lasts until these buffers have been saturated. Therefore, this result suggests a sudden increase in the available buffer concentrations as oscillations start. Figure 2. The model parameters can be fitted to approximate different spikes shapes. A, Experimental Ca 2+ transients measured using the dyes OG and TR. The y axis is the ratio OG/TR after the trace has been detrended with the moving average method. The spikes were aligned to produce the least squared difference over the shown range. The average spike is shown for comparison with the model. B, A simulated Ca 2+ transient with buffer parameters optimized to replicate the shape measured with OG. C, An experimental Ca 2+ transient recorded using Cameleon (YC2.1) as a reporter. The y axis is the ratio of YFP/CFP, after the trace has been detrended with the moving average method. The spikes were aligned to produce the least squared difference over the shown range. The average spike is shown for comparison with the model. D, A simulated Ca 2+ transient with buffer parameters optimized to replicate the shape measured with Cameleon. The y axes of the simulated data represent the Ca 2+ bound to the fluorescent dye/Cameleon complex and thus represent what would be measured in the experiments. The buffer parameter values are given in Table III Furthermore, when the simulations are initiated with one buffer (K d = 5.4 mM) and we add in another buffer (k 2 = 0.006 s 21 and K d = 1 mM) during simulation, then we observe high frequency oscillations, much like the early-stage, high-frequency transients, before a new stable oscillatory pattern with a different frequency is reached (Fig. 5A). If we hypothesize that the perception of Nod factor causes an increase in buffering capacity in the nucleoplasm, then the model predicts that a period of rapid oscillations should occur also if further Nod factor is added during existing oscillations. To test this, we reapplied Nod factor to a cell that was already oscillating (Fig. 5B). As can be seen, the effect predicted by the model could be experimentally reproduced. These observations suggest that Nod factor perception may be accompanied by the production or transport of buffers by the cell to the nucleus.
A further observation in the case of multiple buffers concerns signal initiation and termination. We find that increasing the amount of buffer can terminate oscillations (Fig. 6). This observation suggests that secondary buffers and/or the increase in the concentration of a buffer already present may result in signal discontinuation. Either the gain or loss of buffering capacity from the nuclear environment could lead to termination of the Ca 2+ oscillations. Likewise, such changes in buffer concentrations or kinetics could be the essence of signal initiation. A change in the amount of buffers with low dissociation constants can drive the system to oscillate.

DISCUSSION
We developed ODE models of symbiotic Ca 2+ oscillations in legumes to suggest plausible key components and mechanisms for signal generation. The basic model consists of (1) a Ca 2+ activated ion channel whose function is to hyperpolarize the nuclear membrane (DMI1), (2) a membrane voltage-activated Ca 2+ channel that releases Ca 2+ into the nucleoplasm, and (3) a Ca 2+ pump that pumps Ca 2+ back into the nuclear envelope lumen (MCA8). We have shown that the behavior of these three components is sufficient to produce selfsustaining Ca 2+ oscillations.
The inclusion of buffers in our system has pronounced effects, consistent with previous observations (Marhl et al., 1997;Falcke 2004). As shown here, the changes in buffering capacity offer an attractive explanation for both the initiation and termination of oscillations. Furthermore, the model predicted that increasing the buffer capacity of an oscillating system would introduce an intermediate phase of rapid oscillations, which was verified experimentally. The observation of rapid spiking upon application of Nod factor prior to sustained, lower frequency oscillations can therefore be understood, as suggested by the modeling, as a consequence of higher buffer availability. These high frequency transients have been experimentally observed both when the cell has not previously come in contact with Nod factor as well as when the cell is already spiking due to prior application Figure 3. The plot shows how changing the total buffer concentration can increase or decrease the period of the symbiotic Ca 2+ oscillation. The simulation has relatively large dissociation and association rates of the Ca 2+ buffers, which also slow oscillations, and it includes two buffers: one that was held at a constant concentration and one that was varied. For the buffer parameter values, see Table III. of Nod factor. From this we infer that the introduction of Nod factors to the system also brings with it the release of buffers, possibly the migration of calmodulin to the nucleoplasm from the cytosol, as is observed in animal systems when Ca 2+ levels become elevated (Deisseroth et al., 1998), or other changes in the nucleoplasm that free up more buffers. This result suggests that components beyond those normally associated with the Sym pathway play a vital role in shaping the Ca 2+ signature. Also, the cell-to-cell variability that is often observed may be caused by different concentrations and characteristics of buffers in the different cells.
Differences in Ca 2+ oscillation signatures have been seen not only between genetically identical cells, but also between legume species. The aquatic legume Sesbania rostrata has been shown to have faster oscillations when nodulated at lateral root bases, compared with the oscillations associated with root hair invasion in M. truncatula (Capoen et al., 2009). Such differences could be caused by buffers, but could also be caused by other differences, for example, the pump rate. Our model can produce different spike shapes and periods by varying the pump rate (data not shown); however, only a small range of pump rates gives rise to oscillations. Thus, the pump rate is likely to contribute to the oscillatory behavior, but not be the major determining factor for the observed variation. In this article, we therefore focused on the effects of changing buffer characteristics. Ca 2+ reporters are known to have different binding kinetics, and Ca 2+ buffers are likely to change within cells. By changing the corresponding model parameters, we can account for a range of observations; however, we cannot exclude the, possibly combined, effect of other parameter changes as an alternative explanation.
We hypothesized a purely voltage-controlled Ca 2+ channel. Voltage changes in the nuclear membrane of neural cells have been shown to be involved in periodic release of Ca 2+ from the nuclear envelope (Yamashita, 2011). However, there are other mechanisms compatible with the data, such as the more complex option of voltage and ligand gating or purely ligand gating (Dodd et al., 2010). When more information is available on the identity and regulation of the Ca 2+ channel involved in these symbiotic signals, the proposed model can be modified.
Since we are matching the simulated data to experimental data integrated over one entire cell, this Ca 2+ model does not take into account spatial information.  A possibility for further developments would be to extend this to a spatial model, which can then include, for example, the location of channels and diffusion rates. This would also allow for testing of other regulatory mechanisms of channels that require spatial information, such as Ca 2+ -induced Ca 2+ release, which is a common mechanism in animal cells (Meyer and Stryer, 1991).
Recently, a study described symbiotic Ca 2+ oscillations that occur in the cortical cell layer at later stages of symbiotic entry into the plant root (Sieberer et al., 2011). The authors showed that in both rhizobial and arbuscular mycorrhizal symbioses, the path of progression of the symbiont from the epidermal cell layer down to the cortical cells is preceded by slow Ca 2+ oscillations, and, interestingly, the cells show a higher frequency of oscillations as the microbes start to enter the cell (Sieberer et al., 2011). It is interesting to speculate about the cause of this change in oscillation frequency, and one hypothesis could be a change in buffer characteristics in the cell. This could be triggered by a signal from the microbe, such as Nod factors, which would become available as the microbes pass the outer barriers of the cell.
The model presented here is a first attempt to reproduce the symbiotic Ca 2+ oscillations based on three proteins that are genetically required or predicted to be required for the oscillations. Although our goal has been to provide insights in symbiotic Ca 2+ oscillations, our model is not specific to legumes. Our model builds upon established work in pancreatic beta cells and studies into the effects of buffers in other systems. The proposed working model provides insight into the observed variability of oscillations by studying buffer influences, and it enhances our understanding of how the oscillations are produced.

Plant Material and Ca 2+ Imaging
Two approaches for monitoring Ca 2+ in vivo were used. One method is to inject fluorescent Ca 2+ indicator dyes into root hair cells of wild-type Medicago truncatula (Rudolf et al., 2003), following the protocol described by Capoen et al. (2009) for M. truncatula. The dyes are nonratiometric, using the Ca 2+ -sensitive dye OG 488 BAPTA-1-dextran 10,000 K d and the insensitive TRdextran 10,000 K d (Molecular Probes). The resulting time series is the ratio of the two, OG/TR. The second Ca 2+ measurement method is to transform the plants with the Ca 2+ reporter gene YC2.1, which when expressed indicates Ca 2+ levels by fluorescence resonance energy transfer (Miyawaki et al., 1999). A stable line of M. truncatula R108 YC2.1 was used, and the plant growth and imaging followed the protocol described by Miwa et al. (2006). The resulting time series consist of the ratio of yellow fluorescent protein (YFP) and cyan fluorescent protein (CFP).

Data Processing
Data processing was done in the statistical software R before comparing the experimental data to the simulated data. Large jumps caused by the imaging process were removed. The Ca 2+ time series were detrended using the moving average method (Ghil et al., 2002) to remove the effect of photobleaching. All experimental traces shown are the ratios of OG/TR or YFP/CFP after detrending.

Parameter Estimations
Initial estimates of the parameters were determined such that the concentration profile matched a single Nod factor-induced spike (Supplemental Fig.  S6) using a single shooting algorithm with a particle swarm optimizer with 10,000 particles (He et al., 2007). These parameters are listed and described in Table II. Further parameter refinement was carried out using the following methodology. We extracted all spikes, aligned them so that their squared differences were minimized, and used simulated annealing (Press, 2002) to minimize the least squares distance between the solution to the ODEs and the aligned spikes. We carried out this process for time series obtained using both YC2.1 and OG/TR. These fitted buffer parameters are shown in Table III. The buffer parameters for the supplemental figures are listed in Supplemental  Table S1.

Bifurcation Analysis
The bifurcation diagrams were computed using the software package AUTO in XPP (Ermentrout, 2002). The computations were started from a steady state found by varying the parameter E ps , which relates total Ca 2+ to free Ca 2+ . In some cases, AUTO was unable to trace the limit cycles for the entire chosen region of parameter space because of the computationally intensive calculation required.

Supplemental Data
The following materials are available in the online version of this article.
Supplemental Figure S1. The key symbiotic calcium spiking components and their location.
Supplemental Figure S2. Mimicking the calcium oscillations with buffers.
Supplemental Figure S3. System dynamics of buffers and calcium oscillations.
Supplemental Figure S4. Buffers can give rise to a new limit cycle.
Supplemental Figure S5. Key electrical characteristics as a function of the membrane voltage.
Supplemental Figure S6. Spike shapes can be matched to estimate parameters.