Experimental induction and mathematical modeling of Ca2+ dynamics in rat round spermatids

ABSTRACT Cytosolic Ca2+ concentration ([Ca2+]) has an important role in spermatozoa and hence it regulates fertilization. In male germinal cells, there are indirect evidences that this ion could regulate physiological processes in spermatogenesis. Since little is known about Ca2+ homeostasis in spermatogenic cells, in this work we propose a mathematical model that accounts for experimental [Ca2+] dynamics triggered by blockade of the SERCA transport ATPase with thapsigargin in round rat spermatids, without external Ca2+ and with different extracellular lactate concentrations. The model included three homogeneous calcium compartments and Ca2+-ATPase activities sensitive and insensitive to thapsigargin, and it adjusted satisfactorily the experimental calcium dynamic data. Moreover, an extended version of the model satisfactorily adjusted the stationary states of calcium modulated by extracellular lactate, which is consistent with the participation of a low affinity lactate transporter and further lactate metabolism in these cells. Further studies and modeling would be necessary to shed some light into the relation between Ca2+-lactate-ATP homeostasis and cell–cell interactions in the seminiferous tubules that are expected to modulate Ca2+ dynamics by hormonal factors or energetic substrates in meiotic and postmeiotic spermatogenic cells.

In this work, using the available evidence regarding Ca 2+ homeostatic mechanisms in spermatogenic cells and perturbations in ½Ca 2þ � cyt by the blockade of the SERCA transport ATPase with thapsigargin, we developed mathematical models that describe the dynamics of ½Ca 2þ � cyt in round spermatids without extracellular Ca 2þ . These models are useful in interpreting and predicting ½Ca 2þ � cyt dynamics under physiological conditions or pharmacological interventions.

Animals
Adult (40-60-days old) male Sprague-Dawley rats were acquired from the Animal Facility in the Faculty of Sciences, Universidad de Valparaíso, Chile. The rats were housed in groups of 3-4 animals per cage under a 12 h light:12 h dark cycle with water and rat chow ad libitum. The following procedures were initiated by randomly choosing 1 rat that was euthanized by cervical dislocation following anesthesia with CO 2 . Usually, rats were euthanized in the Animal Facility, and the testes were obtained in a room specially equipped for this procedure.

Ethical statement
All the experiments were conducted in accordance with the rules established by the Consortium for Development of a Guide for the Care and Use of Agricultural Animals in Agricultural Research and Teaching and by the National Research Council. All experimental protocols were reviewed and approved by the Chilean National Fund for Science and Technology (FONDECYT), and the Ethics Committee of the Pontificia Universidad Católica de Valparaíso (EC-PUCV-10/2013). None of the authors have served in this committee.

Rat round spermatid cell isolation
Rat round spermatid cell population was isolated using velocity sedimentation separation in a 2-4% BSA gradient, as described by Romrell et al. [37]. The round spermatid fractions (92 ± 4% purity) were identified both by their size and typical nuclear features after staining with Hoechst 33342 [38]. Most of the measurements in this study were obtained using cell suspensions of round spermatids loaded with specific fluorescent probes, with a Fluoromax 2 fluorimeter (Jobin-Yvon-Spex, NJ, USA).

Intracellular Ca 2þ measurements in rat round spermatids in suspension
Rat round spermatids were suspended (20 × 10 6 cells/mL) in KH medium (140 mM NaCl, 4 mM KCl, 1.6 mM MgCl 2 , 1.6 mM KH 2 PO 4 , 10 mM HEPES), pH 7.4, containing or lacking Ca 2þ (0.5 mM CaCl 2 or 1 mM EGTA, respectively) and having either 5 mM L-Lactate (KH-lactate), or 5 mM glucose (KH-glucose). The cells in KHlactate were loaded with 5 µM of the calcium probe Fura-2 AM by incubation for 1 h at room temperature under an O 2 atmosphere, then washed three times at 4°C, and finally re-suspended in the same medium. Fluorescence measurements were performed after adding a concentrated cell suspension (50 µl) to a temperature-regulated and stirred spectrofluorometer cuvette (3.0 mL) that contained 2.55 mL of the four different media described above, giving cell concentrations of 2.94 × 10 6 cells/ mL, to which thapsigargin was added. The ½Ca 2þ � cyt determinations were performed using a ratiometric method as described by 39, (39]. Fura-2 calibration was performed by cell lysis with digitonin (20--25 µg/mL) and further addition of 1 mM EGTA to determine Fmin, and addition of 3 mM CaCl 2 to the digitonin-treated cell suspension to determine Fmax. All measurements were performed without adding external Ca 2+ and 1 mM EGTA, unless indicated otherwise.

Effects of thapsigargin on lactate-induced ½Ca 2þ � cyt
Cells were treated with two lactate concentrations known to induce different basal ½Ca 2þ � cyt levels [5], for 0, 5, and 10 min, after which thapsigargin (300 nM) was added. The dynamics of ½Ca 2þ � cyt estimated using fura-2 measurements are shown in Figure 1, and these are similar under different conditions. However, there is a slight tendency for ½Ca 2þ � cyt to be higher at longer times of lactate incubation before thapsigargin addition. Figure 2 shows measurements of ½Ca 2þ � cyt for greater range of extracellular lactate concentrations. Because the steady-state condition was not reached in these measurements, the steady-state ½Ca 2þ � cyt value was obtained by extrapolation using exponential decay. Table 1 shows a decreasing relationship between extracellular lactate concentration ([L]) and ½Ca 2þ � cyt before thapsigargin addition and at the extrapolated steady-state ½Ca 2þ � cyt after thapsigargin addition. The data were obtained from triplicate experiments similar to Figure 2.
As shown by Herrera et al. [5], ½Ca 2þ � cyt was dependent on [L] in round spermatids. The dependence of the post-thapsigargin steady-state ½Ca 2þ � cyt value on [L] strongly suggests that homeostatic mechanisms other than SERCA-ATPases operate following lactate (and its metabolism [5]) treatments of these cells. As will be shown further in this article, these results can be used to propose a mathematical model linking energy metabolism and ½Ca 2þ � cyt homeostasis. Mathematical models to describe thapsigargin-induced ½Ca 2þ � cyt dynamics and lactate-modulated steady-states ½Ca 2þ � cyt in round spermatids Figure 3 shows a schematic of a round spermatid composed of three intracellular Ca 2+ distribution compartments: cytoplasm, endoplasmic reticulum (ER), and acidic vesicles (AV). In this Figure the following calcium fluxes are proposed: J 1 : Ca 2þ entry from the cytoplasm to the ER due to a SERCA type Ca 2þ ATPase that transports 2 Ca 2+ ions per hydrolyzed ATP [40][41][42][43].
J 2 : Ca 2þ entry from the cytoplasm to AV due to a Ca 2+ -ATPase that transports one Ca 2þ per ATP hydrolyzed [44].
J 3 : Ca 2þ leakage from the ER lumen toward the cytoplasm.
J 4 : Ca 2þ leakage from the lumen AV toward the cytoplasm. J 5: Ca 2þ efflux from the cytoplasm toward the extracellular medium due to Na þ =Ca 2þ exchanger [45].

Basal mathematical model for thapsigargin-induced ½Ca 2þ � cyt dynamics
We searched for a minimum number of flux components to model the ½Ca 2þ � cyt dynamics in round spermatids after thapsigargin addition, and the following were the preliminary assumptions: 1. The total calcium in the cell was assumed constant (e.g [46][47][48][49][50][51]), that is, the cell behaves kinetically as a closed system for Ca 2+ . This is because outside the cell there is very low free ½Ca 2þ � (estimated [Ca 2+ ] ext = 5 nM) and the Ca 2þ efflux from the cell is assumed to be null (J 5 = 0).
2. In the cytosol, with volume Vol cyt , there is an f cyt ratio of free Ca 2þ to total calcium in the same compartment. Similarly, in the ER, with Vol ret , there is an f ret ratio of free Ca 2þ to total calcium in the ER; and in AV, with volume Vol ves , there is an f ves ratio of free Ca 2þ to total calcium in AV.
3. J 3 and J 4 are proportional to the electrochemical gradient across the membrane.
4. For J 3 , the potential difference is constant [52] and is equal to 0 mV [53]. 5. For J 4 , it is assumed that the membrane potential of the compartment is not zero and meets the Goldman-Hodgkin-Katz equation [54].  6. The membrane area of the different compartments is considered constant during the experiment [55].
According to the previous definitions of fluxes in the scheme of Figure 3, the net fluxes inside the cytosolic, reticular and vesicular compartments (named J Net Cacyt ; J Net Caret and J Net Caves ; respectively) are given by and its relationship with the derivation of the concentration of Ca 2þ in each homogeneous compartment is given by [56] J Net Since we have hypothesized that the cell does not exchange Ca 2þ with the extracellular compartment, the total intracellular Ca 2+ must be conserved as follows: Additionally, in this part, we assume that thapsigargin has been applied, then J 1 ¼ 0. Moreover, from Eqs. 1-7 and A10-A12 (with the mathematical expressions for the fluxes, according to Appendix A), the following relations were obtained Such that

Basal mathematical model was adjusted to the experimental data by applying 0.3 μM of thapsigargin to cells after 5 min of incubation with 10 mM lactate
Before adding thapsigargin, ½Ca 2þ � cyt values in rat round spermatids were recorded in the ranger of 15 to 50 nM. After 5 min of incubation, 300 nM thapsigargin was added, causing the rapid rise of ½Ca 2þ � cyt to a maximum between 300 and 500 nM, and subsequently decreased in slow, exponential, and asymptotic mode, reaching values between 100 and 200 nM. However, these stationary values were not reached during the experiments and, instead, were estimated theoretically by an exponential decay extrapolation. Figure 4 shows the results of 15 different experiments, including all the theoretical fit by means of the differential equation model, as described below.
Using the function ode45 from the MATLAB software, the theoretical values of ½Ca 2þ � cyt were obtained by the fifth order Runge-Kutta method. The values of the adjustment parameters for each experiment ( Table 2) were obtained using MATLAB's function fminsearch, which applies the Nelder-Mead method [57,58], to minimize the mean squared error given by the differences between the theoretical and experimental values ½Ca 2þ � cyt for the total recording after adding thapsigargin. Table 3 shows the statistical analysis of the obtained parameter values, considering dispersions, averages, and the initial parameter values used to start the adjustment. Initial parameters were calculated from the corresponding definition equations using the parameters listed in Table 4.

Modeling lactate-modulated steady-states ½Ca 2þ � cyt
The following is a proposed model for the time change of the cytoplasm [ATP]: where L ½ � is considered constant, ν a constitutive ATP synthesis rate independent of L ½ �; f an ATP synthesis function dependent on L ½ �; and k, a firstorder constant of ATP hydrolysis, a simplified way of representing the influence of several ATPhydrolyzing cell processes. The function f was associated with lactate transport as a limiting factor of the metabolic flux of lactate, and the transport activity was described by a rectangular hyperbola including an affinity constant K L : However, assuming a steady-state for Ca 2þ ½ �cyt.
The ½Ca 2þ � cyt at steady-state (½Ca 2þ � cyt st ) can be calculated as a function of [L] from Eqs. 8-9, 13, 17-19 y A7-A8 Table 2. Estimated parameters obtained from the mathematical model adjustment. The first column corresponds to each experiment shown in Figure 4, giving 15 sets of adjustment parameters.
With thapsigargin. and having defined: In the absence of thapsigargin, besides the AV Ca 2þ -ATPase pump, we have activity of the ER Ca 2þ -ATPase pump (J1 > 0), and similar calculations allow us to obtain an equation such as Eq. 20, but replacing δ with δ 0 : with Without thapsigargin.
Other parameters are defined in Appendix A. In this case, as a state immediately prior to the application of thapsigargin (Figure 2), the assumption of the steady-state may not be adequate, being rather a first approximation. The proposed reason for deviations from noncompliance with the steady-state in the referred case will be discussed later.
For simplicity, the K MATP affinity constants of ATP for ER and AV Ca 2+ -ATPases were assumed Table 3. Analysis of the adjusted parameters shown in Table 2. It is shown the average ( � PÞ, the coefficient of variation as percentage (CV%, equal to 100xSD/ � P) and the ratio between the average values of the adjusted parameter to the initial estimated value of the parameter ( � P/P i ), for each adjusted parameter.  to be the same. However, both the number of each pump and their maximum rates may be different for each case, which is indicated by the primary and non-primary parameters (Appendix A). The model adjustments to the steady-state, given by Eq. 20 and 25, are shown in Figure 6. The adjustment parameters are shown in Table 5. So far, the integrated Ca 2+ flux model represented by the dynamic system given by Eq. 8 and 9, was adjusted to the experimental data shown in Figure 4. In its extended version, that is, by considering the functional dependence of ½Ca 2þ � cyt on extracellular [L], the model was also adjusted to the data shown in Figure 5. Next, we explored the possibility of using alternative models.

Other models for adjusting the dynamic and steady-state calcium data
We tested a simpler model that does not consider Ca 2+ -ATPases in AV. However, in such a model, the stationary states of ½Ca 2þ � cyt in the presence of thapsigargin are independent of ATP, and hence it would be lactate independent, unlike what is observed in the experimental data (Figure 5b) Similarly, the change of the original closed model into an open model using a Na þ =Ca 2þ exchanger in the plasma membrane (Appendix B) could not be adjusted to the data in Figure 4. In this model, the steady-state ½Ca 2þ � cyt after thapsigargin does not depend on [ATP] (Eqs. B21 and B22); therefore, it cannot be adjusted to the data points shown in Figure 5b.

Discussion
In this work, the cytosolic calcium of round rat spermatids was studied by fluorometric methods. Since the measurements were applied to samples composed of many millions of cells, it is possible that the dynamics of calcium per cell are qualitatively different from the average dynamics obtained in these experiments. In any case, the basal calcium records obtained were stationary or slightly increasing. Following the addition of thapsigargin, an abrupt increase and a slow exponential decay toward a steady-state were observed. These records were modeled mathematically considering the structural elements of a likely average model cell. With respect to obtaining parameters to adjust the proposed model to experimental data, it is important to mention that the parameters used at the beginning of all the adjustment algorithms were determined according to bibliographic data, experimental evidence, or Table 5. Adjusted parameters for Equation 20 and 26. The adjusting parameters were obtained from a descending-step algorithm from initial values. Without thapsigargin (Figure 5a), the initial values for the adjusting parameters were δ 0 0 =δ 0 and Φ 0 , fixing the other parameters to those found for the adjustment in Figure 5b. With thapsigargin (Figure 5b), the initial values for the adjusting parameters were δ 0 =1.67 ×10 −4 (d/p), Φ=0.3162, ω 0 =1.4142, K L0 ¼ 2, B 0 ¼ 1.2083 x10 −4 ((q-r)/p).  Table 1. All the adjusted parameters are shown in Table 5.
simple convenience to obtain biologically plausible parameters (Appendix C). The proposed mathematical model to reproduce the experimental dynamics of ½Ca 2þ � cyt in round rat spermatids considers the calcium fluxes given in Figure 4 and structural elements such as the number of compartments, their homogeneity, their volume and their free calcium fractions, in addition to the parameters of each Ca 2þ flux. The first Ca 2þ measurements were made to optimize the recording time and lactate concentrations used. Thapsigargin, a Ca 2þ -ATPase 1 blocker present in the ER, was used to perturb the cytosolic calcium and test the relevance of the mathematical model.
In the experiments where thapsigargin was added to cells incubated with lactate, different records of ½Ca 2þ � cyt versus time were obtained. However, despite being different cell samples, in all data experiment the same set of parameter values was used to initialize the adjustment algorithm. Finally, after parameter adjustment, a different set of optimized parameter values was obtained for each experiment. This diversity of parameter settings is consistent with the fact that each experiment was conducted on different samples (Figure 4, Tables 2 and 3). Table 3 shows that the parameters obtained in Figure 4 have standard percentage deviations between 26% and 107%. However, comparing the averages of each obtained parameter values with respect to the initially assumed parameter values (that is, the values to initialize the adjustment algorithm), it was observed that the values of g were very similar, while p, r, and d did not differ by more than one order of magnitude, and q differed by 5 orders of magnitude. In the latter case, it could be that the complex composition of q by other parameters (Eq. 11) determines a highly amplified global error. Moreover, in Eq. 8 we must consider that q[Ca 2+ ] cyt is much less than d[Ca 2+ ] 2 . Therefore, this difference (corresponding to at least 5 orders of magnitude, according to Table 3 and Figure 4) would determine an increased sensitivity of the q parameter with respect to measurement and rounding errors. On the contrary, the aforementioned similarity between the adjusted g and the assumed initial g would agree with the composition of the parameter g (Eq. 14) being simpler and that the assumed values of the parameters that compose it are more reliable ( f ret ; A ret ; Vol ret ). In any case, the proposed model was adequately adjusted to each of the experimental records shown in Figure 4. Although all adjustment parameter values are within acceptable physiological ranges (Table 3), it is possible that for any of the 15 experimental records in Figure 4 the obtained set of parameter values is not unique. Therefore, the current values should be verified in future experiments, such as for single cell recording experiments. Any suitable mathematical model of the dynamics of ½Ca þ2 � cyt should be able to explain the dependence of stationary calcium states on lactate (Statistical adjustment in Figure 3, adjustment of the model proposed in Figure 6). In that sense, a simpler model, without Ca 2þ -ATPase, does not account for such stationary states with thapsigargin, nor a more complex model, with Na þ =Ca 2þ exchanger activity. Instead, using the model originally proposed in this work, we observed that by replacing some parameters by expressing ATP ½ � as a function of L ½ �, it was possible to model ½Ca 2þ � cyt stationary states before thapsigargin and after its application. It was assumed that the rate of ATP synthesis is given by two terms: a constant term that represents constitutive synthesis of ATP that is independent of lactate (ν) plus one term proportional to a hyperbolic function of extracellular lactate [59]. Furthermore, it was assumed that the overall ATP degradation process is represented as first-order expression. In this extended model, which incorporates L ½ � to explain variations of ½Ca 2þ � cyt stationary states, it was assumed that ATP affinity constants for both Ca 2þ -ATPase (1 and 2, ER and AV, respectively) were equal. As a result of the adjustment of this extended model to the data of stationary states (pre-and post-thapsigargin), it was found that the set of adjustment parameters meets the following conditions that validate the model: -Parameter difference: δ 0 > δ, which indicates that stationary states are more influenced by Ca 2þ -ATPase activity in case without thapsigargin than with thapsigargin. In particular, of the values shown in Table 5, we can assume that Ca 2þ -ATPase activity in the ER membrane is approximately 17 times that of the AV membranes (ðδ 0 À δÞ=δ ¼ 17Þ. -Equal parameters: The following parameters are similar in cases of steady-state pre and postthapsigargin ( Figure 6 and Table 5): B (includes parameters that determine the dynamics of calcium), ω (includes parameters of transport and lactate metabolism) and K L (Affinity constant of the lactate transporter). Interestingly, a K L ¼ 53 mM was obtained, a high Km value as reported for the low affinity lactate transporter [59,60]. Regarding Φ ¼ ν kK MATP , our model requires the same parameter values in the stationary states before and after thapsigargin treatment. However, Φ diminished 24% after thapsigargin addition. This decrease is interesting, because, it could be due to a hypothetical decrease in ν, that is, a diminished constitutive synthesis of ATP. Moreover, this is consistent with the notorious deviation from the steady-state condition prior to thapsigargin in cases of low L ½ �, where there is a very low hyperbolic term contribution to the synthesis of ATP.
In conclusion, in the proposed mathematical model for round rat spermatids, consisting of a closed system, three homogeneous compartments, and Ca 2þ -ATPase activities sensitive and insensitive to thapsigargin, adjust satisfactorily to the experimental calcium dynamics data. In its extended version, which considers elements of energy metabolism, it also adjusts satisfactorily to the stationary states of calcium modulated by lactate, while a spontaneous decrease in constitutive synthesis of ATP would be enough to explain some deviations from the proposed mathematical model. The compartments of the model can be the sum of several compartments of the same type. In another area, the adjustment of the model to the steady-state data modulated by lactate in round rat spermatids suggests the involvement of a lowaffinity lactate transporter, which would be a limiting element in the synthesis of ATP for the activity of the Ca 2þ -ATPase activities AV pumps. Therefore, more research is necessary in this regard. In broader terms, further studies would be useful to clarify important details of lactate transport and ATP homeostasis, so that a more precise mathematical model can be developed to contribute to the study of intracellular calcium modulated by hormonal factors or metabolic energy substrates.