Modular analysis of the control of flagellar Ca2+-spike trains produced by CatSper and CaV channels in sea urchin sperm

Intracellular calcium ([Ca2+]i) is a basic and ubiquitous cellular signal controlling a wide variety of biological processes. A remarkable example is the steering of sea urchin spermatozoa towards the conspecific egg by a spatially and temporally orchestrated series of [Ca2+]i spikes. Although this process has been an experimental paradigm for reproduction and sperm chemotaxis studies, the composition and regulation of the signalling network underlying the cytosolic calcium fluctuations are hitherto not fully understood. Here, we used a differential equations model of the signalling network to assess which set of channels can explain the characteristic envelope and temporal organisation of the [Ca2+]i-spike trains. The signalling network comprises an initial membrane hyperpolarisation produced by an Upstream module triggered by the egg-released chemoattractant peptide, via receptor activation, cGMP synthesis and decay. Followed by downstream modules leading to intraflagellar pH (pHi), voltage and [Ca2+]i fluctuations. The Upstream module outputs were fitted to kinetic data on cGMP activity and early membrane potential changes measured in bulk cell populations. Two candidate modules featuring voltage-dependent Ca2+-channels link these outputs to the downstream dynamics and can independently explain the typical decaying envelope and the progressive spacing of the spikes. In the first module, [Ca2+]i-spike trains require the concerted action of a classical CaV-like channel and a potassium channel, BK (Slo1), whereas the second module relies on pHi-dependent CatSper dynamics articulated with voltage-dependent neutral sodium-proton exchanger (NHE). We analysed the dynamics of these two modules alone and in mixed scenarios. We show that the [Ca2+]i dynamics observed experimentally after sustained alkalinisation can be reproduced by a model featuring the CatSper and NHE module but not by those including the pH-independent CaV and BK module or proportionate mixed scenarios. We conclude in favour of the module containing CatSper and NHE and highlight experimentally testable predictions that would corroborate this conclusion.


Introduction
Intracellular calcium ([Ca 2+ ] i ) plays a key role as a second messenger in a wide variety of biological processes. Its elevation above resting levels conveys diverse physiological signals and may trigger cellular processes such as programmed cell death and egg activation. Both versatility and specificity of this ion as a physiological signal lie in how it is temporally and spatially organised within a cell [1]. Calcium oscillations are widely described in both excitable and non-excitable cells, and lie at the core of several signalling processes, such as neuronal firing, embryonic cell differentiation, immune cell activation and rhythmic beating of the heart. The mechanistic basis of these oscillations have been worked out in a number of cells. A good example is the calcium induced calcium release, which is generated by an inositol trisphosphate and calcium signalling system that controls a host of biological processes [2], for which its operation has been understood with the help of mathematical models [3,4]. Furthermore, in olfactory neurons, calcium influx is regulated directly by cyclic nucleotides [5,6].
The guidance of sea urchin spermatozoa towards their conspecific eggs during broadcast spawning events is orchestrated by [Ca 2+ ] i -spike trains elicited by the chemoattractant sperm activating peptides (SAP) released by the eggs. These spike trains are decoded into a series of acute turns followed by straight swimming episodes that define the trajectory of the sperm cell. Preventing calcium influx from the external medium will abrogate the cell's capacity to orientate in space. The sequence of spikes can be triggered by SAP in sperm that are fixed, showing a typical temporal organisation, characterised by a gradual decrease of the amplitude of the spikes and progressive increase of the interspike intervals.
The signalling network that mediates calcium-steered chemotaxis in sea urchin sperm has been an object of intense study, and several decades of experimental research have implicated multiple channels and molecular components [7,8]. It is now clear that the first step leading to Ca 2+ influx is an increase of intracellular cGMP concentration, which is synthesised in response to the natural ligand. Experiments in which the receptor activation was bypassed by releasing a caged analogue demonstrated that cGMP is sufficient to trigger the signalling network downstream [9]. The rise in cGMP leads to membrane hyperpolarisation by opening a cGMP-gated K + -channel (KCNG). How the hyperpolarisation is then transduced into Ca 2+ flux dynamics is not fully established. It has been hypothesised [10,11] that the membrane potential shift to more negative values would have two effects: to stimulate a depolarising cationic current by sperm specific hyperpolarisation-activated cyclic nucleotide-gated channels (spHCN) and to remove the inactivation of voltage-gated calcium channels (Ca V ), which would then further depolarise the membrane. The depolarisation of the membrane would close spHCN and Ca V , allowing hyperpolarisation fluxes to predominate, repeating the cycle. This cycle would be the basis of the [Ca 2+ ] i -spike train, the temporal organisation of which could be modulated by Ca 2+ -activated Cland K + channels. Another hypothesis [12] is that the cGMP-driven hyperpolarisation, by activating a Na + /H + exchanger that extrudes protons, would rise the intraflagellar pH (pH i ) eliciting Ca 2+ influx by a pH-dependent Ca 2+ channel, such as CatSper. This is a sperm specific calcium channel with unique features that separate it from the classical Ca V family, including modulation of voltage-gating by pH i and Ca 2+ . These features open the possibility for another negative feedback underlying [Ca 2+ ] i oscillations. The pH i rise would promote the opening of CatSper that drives an upsurge in [Ca 2+ ] i and membrane depolarisation, which in consequence would inactivate respectively CatSper and the exchanger. The cycle would reinitiate once the membrane hyperpolarises again.
The evidence for the two types of calcium channels proposed to participate in SAP-activated [Ca 2+ ] i -spike trains is summarised in Table 1. Most evidence is indirect and based on pharmacological approaches, which attempt to implicate channels using blocking and activating drugs that often are not specific and therefore have off-target effects. Genetic engineering techniques that would allow definitive conclusions on the necessity and sufficiency of the posited channels (e.g. null mutants) are not available in sea urchin. Therefore, it remains elusive whether a single or multiple types of calcium channels play an active role in sea urchin sperm signalling.
Mathematical modelling may help to explore alternative hypotheses that cannot be directly addressed and disentangled experimentally. The signalling network elicited by Speract in sea urchin sperm has been modelled using the discrete formalism [17,20,18]. While the initial discrete network models featured only Ca V s, a subsequent attempt to gain insights into the effect of multi-target niflumic acid suggested that pH-dependent CatSper would play a key role in Ca 2+ oscillations, and that the calcium-dependent potassium channel (BK) would further modulate the spike trains [20]. The inherent limitations of discrete models in describing concentrations and kinetic details hinder the conclusions that can be drawn from comparison Table 1. Summary of experimental evidence on the main calcium channels proposed in the pathway activated by SAPs. a It is worthwhile to mention that out of the pharmacological blockers listed, none of them are 100% specific to the associated channels. b Antibody raised against a polypeptide of the α subunit in rat, which also recognises the homologous version in mouse and human sperm. c Antibody raised against pore-forming subunit polypeptides of A. punctulata [12] and S. purpuratus [18]. d It is of note that direct electrophysiological recordings for this ion channel have not been achieved in mature sea urchin spermatozoa compared with the homologous complex in mouse, human and macaque spermatozoa [19].
with quantitative experiments. A continuous model was used to identify the core components of the signalling network [21], implicating Ca V channels, but the model was not directly compared to experimental data.
In the present work, we explore which putative calcium channels may mediate the Ca 2+ response to SAP by describing the network using coupled ordinary differential equations that can fully harness the quantitative and kinetic information available in experimental data. We asked which channels may intervene, and how they may control the relevant dynamical features of the observed [Ca 2+ ] i -spike trains. Our aim is to contrast the different hypotheses by their quantitative implications in order to identify critical properties that can be used to experimentally disentangle them.

A modular signalling network
The potential signalling network illustrated in Fig 1 is organised in three modules: an Upstream module, which includes the ligation of the receptor by SAP and the early transient drop of membrane potential V generated by the cGMP response, and two alternative downstream modules that can generate Ca 2+ oscillatory responses. The Ca V +BK module includes voltage-dependent Ca 2+ -channels and Ca 2+ -dependent K + -channel BK; while the CatSper +NHE module is composed of the voltage-, pH i -and Ca 2+ -dependent channel CatSper, the voltage-dependent sodium-proton exchanger and the intraflagellar proton concentration. The variables, parameters and differential equations used to model the components of this network are described in Sec. 4.1 and Sec. 4.2.
The modular structure facilitated modelling and analysis. As a strategic simplification, we hypothesised that most of the components of the Upstream module were independent of the variables downstream and studied it isolated from the remaining dynamics. We tested the validity of this simplifying assumption by showing that this Upstream module accounted quantitatively for experimental time series of cGMP and membrane potential in response to SAP (Sec. 2.2). Once the values of the parameters of this module were obtained, they remained fixed in the subsequent analyses of the dynamics of the downstream components. In these analyses, the modules featuring Ca 2+ -channels were constrained such that given the input of the Upstream module, they generated an output that was consistent with the observed Ca 2+spike trains, i.e. these two putative modules were constrained by independent data sets at the both input and output levels. Hence, by coupling the downstream Ca V +BK and CatSper+NHE modules to this Upstream module, we show in Sec. 2.3 that either module can independently account for the observed Ca 2+ -spike trains patterns.

The Upstream module of the SAP-activated response explains the long term kinetics of cGMP independently of the downstream signals
cGMP is the core component of the Upstream module (Fig 1). cGMP levels rise in response to SAP-induced guanylate cyclase activity and activates KCNG channels that in turn lead to hyperpolarising currents. A. punctulata sperm stimulated with the SAP Resact exhibit a rapid cGMP response that peaks before *400 ms and rapidly drops to a pseudo plateau that slowly decays over several seconds, as observed by stopped-flow kinetics [22]. This biphasic dynamics prompted us to postulate that SAP receptors transit irreversibly through three forms with different levels of associated guanylate cyclase activity, and that the guanylate cyclase and phosphodiesterases activities are insensitive to downstream signals (Fig 1; Sec. 4.2.1 and Sec. 4.1.2). We are thus neglecting putative regulation of phosphodiasterases activity by cGMP itself and by other components of the signalling network, such as pH i or [Ca 2+ ] i [23]. In order to assess these assumptions we fitted a model of the Upstream module to cGMP measurements in bulk cell assays wherein sperm were stimulated with pulses of uncaged SAP along a concentration range spanning several orders of magnitude [22]. The good quality of the fitting (S1 Fig), involving appropriate scaling factors, supported quantitatively the assumptions and provided estimates of kinetic rate constants for receptor state transition, cGMP synthesis and decay (kinetics illustrated in Fig 2A for a SAP concentration that saturates the response of the membrane potential V [10]). The long term dynamics of the average cGMP level in the sperm population is described well by our single cell model, which as a first approximation neglected any potential feedbacks from downstream processes. This is our first modelling result.
By adjusting the density of KCNG and spHCN channels to membrane potential from sperm population data [10], we found values that reproduce the early transient drop in The modular organisation of the signalling network transducing SAP signals to Ca 2+ -spike trains. The structure is separated into 3 modules coupled by the membrane potential variable (V). In the Upstream module (described by the variables S, hc and f o hc ) free receptors R F bind SAP molecules and transit irreversibly through three receptor forms, each with less guanylate cyclase activity: High R H , low R L and inactive R I ; the cGMP synthesised by these receptors binds and opens KCNG channels, which conduct a hyperpolarising potassium current. The membrane hyperpolarisation promotes the opening of spHCN channels, which exert the opposite action on V when conducting a cationic inward current. Two candidate modules are presented to explain the calcium spike trains: one that includes classic voltage-dependent Ca 2+ channels and BK channels (described by variables f c cv , f o cv , and f o bk ), and another that considers CatSper, NHE exchanger and proton concentration (described by f o cs ¼ m cs h cs , f a nh and H). Note in the Ca V +BK module, that the calcium channel is purely voltage dependent, whereas CatSper in the other module has threefold regulation, namely Ca 2+ , H + and V. The complexity of the CatSper channel is illustrated with two gates. Voltage-regulated channel transitions are depicted with blue arrows, where the darker or lighter tones indicate whether the transition is promoted by membrane hyperpolarisation or depolarisation respectively. The blue dotted lines linking some of these arrows to the membrane potential box top or bottom provide the same information. The components in the middle dashed boxes, namely cGMP (G), Ca 2+ (C), proton (H) and membrane potential (V) are the main experimental observables. The ions K + and Na + are depicted in grey to indicate that the model describes only the currents while neglecting concentration changes. membrane potential and the early recovery (Fig 2). To illustrate the separate contribution of each channel, we showed a reference case where KCNG is the only channel present (dotted line in Fig 2) and, as expected, its opening alone drives membrane potential towards E K � −80 mV (potassium Nernst potential). The latter reference value sets a lower limit to the V response that is well below the value reached at the minimum V produced by the KCNG +spHCN case. It is worth noticing that under this scenario, in which only these two channels are assumed to be present, it is not possible to find a given combination of channel densities such that the spHCN current alone counterbalances the hyperpolarising current carried by KCNG, so as to bring membrane potential back to resting levels without abrogating the hyperpolarisation pulse. Therefore, full recovery of voltage requires the additional activity of other ion channels with depolarising currents after these initial events (namely the calcium channels).

Two alternative signalling modules can explain semiquantitatively the structure of the Ca 2+ -spike trains
The [Ca 2+ ] i spike trains observed in individual sperm cells responding to the SAP stimulus are noisy but show some conspicuous trends, as illustrated by those of S. purpuratus sperm ( Fig  3A). The regression coefficients of the relative spike amplitude b A and the relative interspike interval b T over the spike index (illustrated on Fig 3A right) provide numeric measurements of these trends in the individual responses (Sec. 4.4). The amplitude of [Ca 2+ ]-spikes declines progressively in 87% of the responses analysed (Fig 3B top), whereas the interspike interval tends to increase in 65% of the cases (Fig 3B bottom). The bivariate plot of b A vs b T (Fig 3C) shows that a majority of 56% of the [Ca 2+ ] i responses to SAP, as well as the median coefficients, fall within the bottom right quadrant. The convex hull defined by the experimental values that fall within this quadrant (the dashed polygon in Fig 3C) offers a first selection criterium to assess how close the dynamics predicted by the model are to the observations. The complementary criteria are the mean interspike interval and the steady state [Ca 2+ ] i value before the stimulus (Sec. 4.4). The rationale for the first selection criterium is twofold. First, it is the most representative pattern of the experimental spike trains analysed here as well as of those reported in other species (e.g. see fig. 3 in reference [9] for A. punctulata sperm). Second, it is the trend that was most difficult to obtain in the model under either the Ca V +BK or the Cat-Sper+NHE modules, therefore being most suitable to challenge network composition. The spike trains falling in the bottom left quadrant, besides being less frequent experimentally, were easily obtained with the model and therefore were less challenging.
We found that the two alternative signalling modules, Ca V +BK and CatSper+NHE (Fig 1), can explain the characteristic envelope and increasing delay of the Ca 2+ -spikes when coupled alone to the Upstream module. These modules are composed of specific combinations of ion channels and transporters. Each of these combinations has the capacity to produce, given the input of the cGMP kinetics used to calibrate the Upstream module (Fig 2), a series of Ca 2+spikes that are coherent with the predominantly observed ones (Sec. 4.4). The corresponding numerical solutions (Fig 3D and 3E top) are characterised by amplitude and interspike interval regression coefficients (Fig 3D and 3E bottom) that fall within the values of the most likely observed responses in individual sperm cells of S. purpuratus (represented by the labelled circles in Fig 3C).
The two modules and their properties are analysed separately in the following sections. For each module, a first section characterises the Ca 2+ oscillations and explains how they are generated (Sec. 2.4, Sec. 2.6), and a second focuses on explaining the progressive amplitude decay and delay between spikes (Sec. 2.5, Sec. 2.7). We will build on this knowledge to analyse and sperm bound to a coverslip following uncaging of Speract analogue by UV light. The raw measurements (grey) were smoothed using cubic spline (black). The Ca-spikes (vertical grey dashed lines) were calculated as the maxima in the smooth curve higher than 3x the standard deviation of the measurements preuncaging (horizontal gray dashed line). For each spike indexed k, the amplitude A k and the interspike interval T k = t k+1 − t k were calculated. The shaded region indicates the UV light pulse. On the right, the relative spike amplitude A k /A 1 (left axis; black) and relative interspike interval T k /T 1 (right axis; grey) are plotted as a function the spike index k for the two examples. The lines are linear regressions of the relative spike amplitude (black) and relative interspike interval (grey) over the spike index k from which the respective regression coefficients b A and b T were obtained. understand mixed scenarios (Sec. 2.8). The Ca V channel was the first Ca 2+ channel to be implicated in SAP signalling in sea urchin sperm literature and following this chronology we will start the analysis of its module.

Characteristics and mechanism of the [Ca 2+ ] i oscillations generated by the Ca V +BK module
The module features the Ca V channels and the Ca 2+ -dependent BK potassium channels ( Fig  1). The Ca V channels transit irreversibly through three forms-inactive, closed and opendefining an ordered cycle, in which the transition rates are controlled by voltage. Only in the open form Ca V s allow inward Ca 2+ currents that tend to depolarise the membrane. Provided that the transition rates between these forms are of similar order of magnitude, these channels alone give rise to a limit cycle, which as we will see has a constant period that is unrealistically short (Sec. 2.5). This disagreement with the observations called into action the other channel in the module, BK, which once coupled with the Ca V results in larger progressively increasing intervals between the peaks and also increases the peak amplitude.
The numerical solution of the variables of the model including the Ca V +BK module under the reference parameters is shown in Fig 4. The uprise and slow decay of cGMP produced by the Upstream module leads to Ca 2+ -spike trains with an envelope that follows closely the cGMP values (top) and increasing interspike intervals. The initial transient hyperpolarisation (bottom) driven by the opening and currents of KCNG (f o kn and I kn ) precedes and triggers the first peak of [Ca 2+ ] i (C). The nadirs of the [Ca 2+ ] i oscillations are well above the basal concentration of this cation and this minimal plateau of about 500 nM is maintained until the Ca 2+ activity vanishes at about 17 s. A similar oscillatory behaviour with similar period but different phases is observed on the fractions of open channels (second row) and the respective currents (third row), as well as the membrane potential (bottom).
The temporal dynamics of the [Ca 2+ ] i fluctuations predicted by this module can be easily understood qualitatively. The hyperpolarisation induced by KCNG enables the transition from inactive to closed forms of the Ca V and the spHCN-mediated increase in membrane potential opens these closed channels. The open Ca V channels will overcome the hyperpolarising current driven by KCNG and thus depolarise the membrane as Ca 2+ ions flow in. The depolarisation of the membrane leads to the inactivation of the Ca V channels with the hyperpolarisation currents predominating again. At the hyperpolarised membrane potentials the inactive channels will again transit to closed form, which eventually will open as described above, reinitiating the cycle. Intraflagellar calcium, by activating the BK channels, adds a strong hyperpolarising current to that of KCNG. This couples the net hyperpolarising currents to the amplitude of the previous Ca 2+ -spike, in such a way the higher the [Ca 2+ ] i peak, the higher the BK-driven hyperpolarising current will be.  (Table 2). On top, the curves represent the numerical solutions of the variable C as a function of time. The Ca-spikes are identified (vertical dashed lines) and indexed by k (labels). The relative spike amplitude and relative interspike interval are plotted as function of spike index k on the bottom graph together with the regression lines used to obtain the regression coefficients b A and b T for the model, which were represented in C. https://doi.org/10.1371/journal.pcbi.1007605.g003

Control of [Ca 2+ ] i -spike amplitude and interspike interval in the Ca V +BK module
The envelope of the [Ca 2+ ] i spikes follows closely that of the cGMP levels upstream. This led us to speculate that the time course of cGMP might offer a way to understand the fluctuations  Table 2  envelope. To gain a better quantitative insight, a variant of the model was made wherein cGMP concentration, G, instead of a variable became a constant parameter, thus ignoring the upstream variables and reducing the state space (Sec. 4.3). In this simplified model, the qualitative changes of the [Ca 2 + ] i and its period at the stationary states were analysed as a function of the G parameter. This is known as bifurcation analysis in dynamical systems theory, which is useful for characterising the global properties of a system of differential equations, in particular the dependence of asymptotic behaviours on a given parameter of interest (bifurcation parameter). These behaviours can correspond to attractors of fixed point (invariant in time) or limit cycles (stationary periodic), which in turn can be stable or unstable, among others. Wherever there is a qualitative alteration of these behaviours by changing the parameter of interest, it is said that a bifurcation has occurred.
When cGMP concentration is varied in this subsystem, the bifurcation diagram shows 5 bifurcation points (Fig 5A, labeled from I to V). In I, by means of a Hopf-type bifurcation, as cGMP concentration increases above the basal level, the system moves from a [Ca 2+ ] i stable stationary state to a limit cycle. When it reaches II, the limit cycle disappears giving rise again to a stable stationary equilibrium by means of an inverse Hopf-type bifurcation. Further on in III, the system undergoes an inverse saddle-node bifurcation, defined over cGMP  (Table 2). Top: the graphs show the C variable ([Ca 2+ ] i ) as a function of the input parameter G (cGMP concentration). The thick continuous lines are stable equilibria, the dashed ones indicate unstable equilibria, the thin continuous lines are the maxima and minima of stable limit cycles. The bifurcation points are marked with Roman numerals (I to V). The numbered grey circles indicate the peak value for consecutive spikes obtained by numerically solving the complete model (see Fig 4). Bottom: the solid lines represent the period of stable limit cycles as a function of the constant input G. The numbered gray circles correspond to the interspike interval in the numerical solutions, as in the top graph. B) Bifurcation diagrams parameterised by maximal BK conductance density (g bk ). Top and bottom: the lines are as in the graphs in (A); for clarity, the minima and maxima of the limit cycles are omitted. The numbers in the range of 0.0 to 1.1 are the fold change factors multiplying the reference value of g bk (e.g. the curves labelled 1.0 coincide with those of (A), whereas the curves marked 0.0 correspond to a cell without BK channels).
https://doi.org/10.1371/journal.pcbi.1007605.g005 concentration increments, in which a stable equilibrium state coalesces with an unstable one. In this way, in traversing the branch over the stable equilibrium state, the system goes to an unstable state. Now, when advancing on this last unstable equilibrium, this turns to a stable one by means of the direct saddle-node bifurcation (IV). The presence of points III and IV suggest the possible existence of a cusp bifurcation defined in a space of higher dimension, which is corroborated in panel B of the same figure. Finally, the last stable equilibrium state, again by means of a direct Hopf-type bifurcation, gives rise to another limit cycle.
The period of the limit cycle in the upper branch increases as the cGMP concentration decreases (Fig 5A, bottom), but this is not the case for the limit cycle in the lower branch, in which the period shows an opposite yet less pronounced trend. In addition, there is a considerable difference between the two limit cycles in terms of their amplitude, demarcated by the minima and maxima of calcium levels in the bifurcation diagram.
To better understand the role of the BK channel on calcium dynamics in a Ca V -dependent scenario, bifurcation diagrams were calculated within the cGMP concentration range considered in Fig 5A, for different density values of BK, represented by g bk (Fig 5B). In B, unlike the upper panel of the Fig 5A, the lower and upper lines representing the limit cycles were omitted for clarity; their respective unstable steady state values are shown instead (dashed lines). In the bifurcation analysis, the extreme case where g bk = 0, that is, the variant of the model that only presents Ca V , predicts the existence of a single branch for the calcium concentration with a limit cycle whose envelope is controlled by cGMP levels but has a constant period of approximately 0.2 s (see the lines labeled 0.0, in both panels of Fig 5). The gradual addition of BK reveals how bistability is generated (Fig 5B), and exhibits the appearance of an upper limit cycle within the range of cGMP concentration values observed in Fig 5A. The behaviour of the periods associated with the different values of g bk is shown in the lower panel of the Fig 5B. In a three-dimensional representation, with g bk as one of the axes, the bifurcation pattern is a cusp with mounted limit cycles. The presence of the new upper periodic attractor allows the nadir of the [Ca 2+ ] i fluctuations to be placed on a plateau above the basal level, and both their amplitude and increasing interspike interval to approach the physiological response. This implies that the aforementioned properties arise from the coupling of the Ca V and BK channels.
The parameters have been chosen such that as the cGMP concentration values transiently exceed the bifurcation point, the system is forced to oscillate in the upper limit cycle and fall into the lower limit cycle after the decrease in cGMP causes the system to pass through the saddle-node bifurcation (IV). This explains the abrupt [Ca 2+ ] i drop and the low amplitude oscillations that persist afterwards, which are particularly visible in the numerical solution of V (Fig  4, bottom). If the value of g bk were increased above the reference value used in panel A (e.g. by multiplying it by a factor of 1.1), the saddle-node bifurcation would occur at a lower cGMP concentration value and these low amplitude fluctuations would not be revealed (see Fig 5B).
Considering the bifurcations of the model obtained with cGMP as control constant, the behaviour of the full model is easy to understand: as SAP binds to its receptor, a sustained yet decreasing cGMP activity leads the system to approach asymptotically the attractors of the model with constant cGMP concentration, such that the system displays fluctuations whose spike amplitude and interspike interval tend to the amplitude and the period of the attractors (compare the numbered gray circles corresponding to the solution of the full system, which takes into account the temporal evolution of cGMP concentration as in Fig 4, to the asymptotic limit cycles for different fixed values in Fig 5A).
It is important to note that at intermediate levels of cGMP, there is a range of values where limit cycles coexist. From what has been reported in the literature, only the dynamics compatible with the upper limit cycle has been experimentally observed. Presumably, even if the small oscillations of the lower limit cycle were present, they would be masked by measurement errors.

Characteristics and mechanism of the [Ca 2+ ] i oscillations generated by the CatSper+NHE module
We now consider the alternative module whose central component is CatSper, the Ca 2+ channel that has recently been implicated in SAP signalling in sea urchin sperm. We assumed here, that CatSper opening relies on two independent gates: one modulated by both pH i and voltage, and the other one responsive to intracellular calcium. The fraction of open channels is given by the product of the fractions of these two open gates (Eq 44). The module includes also as variables the fraction of active forms of the voltage-sensitive Na + /H + exchangers, f a nh , and the intraflagellar concentration of protons, H. The numerical solutions obtained under the same conditions as in the alternative model are shown in Fig 6. The rise and fall of the cGMP produced by the Upstream module leads to a Ca 2+ -spike train with a slow decaying amplitude and an increasing interspike interval (Fig 6, top). The nadir of the Ca 2+ -oscillations is close or even below the resting level, in contrast with the oscillations predicted by the alternative model (Fig 4). A similar oscillatory behaviour with similar period but different phases is observed on the fractions of open channels (Fig 6, second row), the respective currents (third row), as well as the membrane potential and proton concentration (bottom).
Compared to the Ca V +BK module, the pattern of the Ca 2+ oscillations is differently shaped by the CatSper+NHE module. The individual Ca 2+ -spikes are more prolonged and the corresponding nadir is relatively shorter. The decay of the spike train envelope is less pronounced for the reference parameters in the CatSper+NHE module than in the Ca V +BK module, while the increase in the interspike interval is more marked in the former than in the latter. It is remarkable that the CatSper+NHE module does not predict the sustained plateau of the Ca 2+ -oscillation nadirs present in the Ca V +BK module. Finally, another interesting difference between the two modules is the envelope of the voltage spikes whose amplitudes decrease progressively in the Ca V +BK module and increase in the CatSper+NHE module (bottom in Figs 4 and 6).
The basis for the limit cycle in CatSper+NHE module is qualitatively straightforward. The hyperpolarisation increases the activity of the Na + /H + exchanger, this triggers the extrusion of protons raising pH i . As a consequence of this alkalinisation, the voltage sensitivity of CatSper shifts to lower values (see S3D Fig) and some of these channels open. The inward Ca 2+ currents raise the membrane potential which in turn tends to increase the fraction of open channels, in a positive feedback loop. However, the progressive depolarisation of the membrane inhibits the exchanger and the pH i tends to return to its basal level, shifting the CatSper voltage sensitivity back to higher values. This allows the hyperpolarising current to surpass the depolarising current of CatSper (see Fig 6, third row). The cycle restarts when the membrane becomes again hyperpolarised by the still ongoing KCNG current (Fig 6, second and third row), leading to a new round of enhanced exchanger activity, a transient rise in pH i , and the recovery of Cat-Sper currents that eventually overcome the hyperpolarising currents. Although the feedback inhibition of CatSper by Ca 2+ -inactivation could have a potential role in the limit cycle, its contribution is actually negligible under the reference parameter regime.

Control of [Ca 2+ ] i -spike amplitude and interspike interval in the CatSper+NHE module
As in the case of the Ca V +BK module, the bifurcation analysis of the model featuring the Cat-Sper+NHE module with constant input cGMP reveals a rather complex dynamical structure.
The system displays a cusp involving two saddle-node bifurcations, which, in contrast to the module Ca V + BK, is in lower cGMP concentration ranges, such that two [Ca 2+ ] i stable equilibria coexist close to the resting state. At higher cGMP values, there is a regime where the stable equilibrium solution coexists with a stable limit cycle, whose maxima and minima are indicated by a succession of black circles in the top graph in Fig 7. It is worth noticing that the  Table 2.
https://doi.org/10.1371/journal.pcbi.1007605.g006 coexistence of stable attractors within a given range of cGMP concentrations suggests the presence of bifurcation points that determine the transition of one behaviour to the other, at the boundaries of that range. However, for this model, with the available bifurcation analysis software (e.g. XPPAUT [24]), it was not possible to discover these points or the limit cycle. In fact, https://doi.org/10.1371/journal.pcbi.1007605.g007 the stable limit cycle was detected by solving the model numerically, exploring the phase space by randomisation of initial conditions of the system and allowing it to evolve for long times (>5 × 10 5 iterations). Also an exploration of initial conditions between the two attractors, suggests that a possible presence of highly interlaced attraction basins, which makes numerical calculations difficult.
For this scenario, as was the case with the Ca V +BK module, the amplitude decreases and the period of the limit cycle increases (Fig 7, bottom) as the constant input of cGMP returns to basal levels after having reached high values.
The patterns of the [Ca 2+ ] i fluctuations obtained by solving the full model with the CatSper +NHE module (Fig 6) are now easy to understand. As the cGMP concentration increases rapidly from the resting state, the system falls within the basin of attraction of the limit cycle determined for constant G values; and as the cGMP returns slowly to its basal level, the time required to approach those limit cycles decreases. This can be seen in Fig 7 by means of the empty circles superimposed on the bifurcation diagram, which show how the system tends towards this limit cycle with a spike amplitude reduction and an increase in the interspike intervals. The differences between the asymptotic values of amplitude and periods predicted by this variant of the model with fixed values of G (black circles), and those of the actual peaks obtained when G varies in time (numbered empty circles) are more marked in this module than in the case of the Ca V +BK module.

The decaying envelope and temporal organisation of the Ca-spike trains can be obtained in mixed scenarios involving components of the Ca V +BK and CatSper+NHE modules
In the previous sections, we have shown that both alternative modules, once coupled with the Upstream module, can account for the observed data patterns. These two modules share some features and their components could potentially coexist in the flagellum being coupled through the membrane potential, pH i , and [Ca 2+ ] i . As a first approach to investigate how the coupling of the modules would affect the dynamics, we considered a scenario in which the actual composition of the system would be a weighted combination of the two modules by means of weight parameters given by θ cs = θ and θ cs = (1 − θ) (implementation details in Sec. 4.5). For θ = 0, one recovers the submodel with CatSper module alone, for θ = 1 one recovers the submodel featuring Ca V +BK module only, and when θ = 0.5 the total conductance of the components in each module are half of their respective reference values. Taking the CatSper+NHE module we titrated in the components of the Ca V +BK module under several combination schemes. First, we explored the addition of Ca V and BK simultaneously maintaining their conductance ratio as in the original Ca V +BK module (S4 Fig), then we added either channel alone (S5 Fig). This simple analysis showed that the modules cannot be combined freely without severely compromising the dynamics. Only in situations where one of the original modules largely predominates over the other (either θ *1 or θ *0) the mixed models predict the characteristic decay of the amplitude and the increase interspike interval. BK alone added to the CatSper+NHE module will cancel the oscillations at values as small as 10% of its reference value. At even lower densities it is functionally negligible, modulating the average interspike interval or the duration of the whole spike train without modifying the characteristics of the oscillations. These results are not entirely unexpected considering that in each oscillation of the Ca V +BK module the system spends most of the time hyperpolarised, while in the CatSper +NHE it spends a significant fraction of the cycle depolarised, which is an absolute prerequisite for pH i to return to lower levels closing CatSper gate.
These considerations notwithstanding, the initial analysis of the mixed scenarios was done under a single parameter regime in all the remaining parameters, and it is possible that this may have constrained results. To investigate this possibility we explored the parameter space in search for mixed scenarios that would be coherent with the experimental Ca-spike trains based on the criteria defined in Sec. 4.4. The purpose was not to identify which parameters have the greatest impact on the models' quantitative output, for which we could have resorted to Latin hypercube sampling and quantitative sensitivity analysis as done by others (e.g. [25,26]). Our more immediate goal was to relax the parameter settings in order to detect potential mixed scenarios. Because our assessment of the model output is essentially a qualitative (binary) criterium of coherence with observations, the quantitative sensitivity analysis via partial rank correlation coefficients (as in [25]) or other methods (as in [26]) was not feasible. Based on these considerations, we devised a simple algorithm inspired on adaptive evolutionary dynamics to explore the parameter space. This algorithm (Sec. 4.5 and S1 Algorithm) evolved a seeding population of parameter vectors in a regular grid first neutrally and subsequently under stringent binary selection criteria. Variants of the parameter vectors were generated by a stochastic process, based on a random walk in logarithmic parameter space. The initial grid ensured that the abundances of the components of the Ca V +BK and CatSper+NHE modules were systematically explored in broad ranges, between negligible and excessive conductances. The generation of parameter vector variants ensured a random exploration of the parameter space centred on the realistic estimates of the other parameters, to avoid diversions into spurious parameter regimes. The final ensemble of parameter vectors produced by the algorithm contains both a population that was selected for being coherent with experimental data and an unselected population derived from and closely located to the former (Sec. 4.5 and S1 Algorithm). The projection of the parameter vectors in these two populations on the plane (g cv , g cs ) shows 6 clusters, labelled a-f, depicted in Fig 8A. It is rather conspicuous that the cluster including the reference parameter set for the CatSper+NHE module (cluster b) is narrower than the other clusters. The cluster seems to be more sensitive to parameter changes as the within-cluster fraction of selected parameters is 8%. which is about half the average fraction of selected parameter vectors in the ensemble of the other clusters (16%). The mapping to the plane of the predicted spike amplitude and interspike interval regression coefficients indicates a narrow range of spike amplitude decay ( Fig 8B). The other five clusters, including all mixed scenarios and the cluster centred on the reference Ca V +BK module, are more spread both in parameter space (g cv , g cs ) and in the plane (b T , b A ) summarising the predicted output. To understand the mixed scenarios we asked whether the presence of all the channels Ca V , BK or CatSper was essential to produce the predicted dynamics. To this end, we compared the numerical solution obtained with a given parameter set with the solutions obtained by specifically cancelling the conductance of Ca V and BK on the one hand, or that of CatSper, on the other hand (examples of these numerical solutions are shown in Fig 8 for parameters singled out from each cluster). These comparisons revealed that CatSper has no functional consequences in the parameter vectors in clusters a, c and f whereas Ca V and BK are irrelevant in those of cluster b. In clusters d and e, which include conductance densities of CatSper higher than the reference ones, the characteristic dynamics were severely distorted if either CatSper or Ca V and BK were eliminated. In other words, the reference parameter set for the Ca V +BK module captures the essential qualitative features of the dynamics associated with clusters a, c and f and the reference CatSper+BK module represents well those of cluster b. However, the dynamics corresponding to clusters d and e may not be fully captured by the reference modules.
It is worth noticing that the numerical solutions obtained with the parameter vectors in the CatSper+NHE module cluster (cluster b) show in each oscillation cycle the typical preponderance of depolarisation, with blunt Ca-spikes and marked pH i oscillations. The parameter vectors in the other clusters display very transient depolarisation times, with sharp Ca-spikes ( Fig 8B) and pH i oscillations with shallow amplitude around a relatively high value, qualitatively similar to those obtained with the reference Ca V +BK module (Fig 4). As we will A large number of parameter vectors resulting from a search of parameter space for coherence with experimental time series is illustrated in this figure according to clustering in the projection plane (g cv , g cs ), predicted dynamics, and predicted functionality of the components of the two modules. During the exploration, coherence with experimental time series required that the resting C(0) is about 100nM, the average spike interval is in the range 0.37 s-1.38 s, the predicted (b T , b A ) point falls within the convex hull of most representative experimental values in the bottom right quadrant (polygon in A, right). A) On the left, the clustering of the parameter vectors in projection in the (g cv , g cs ) plane is shown. On the right, the envelope and temporal organisation of Ca-spike trains predicted by each of these parameter vectors is represented as the bivariate plot of the spike-amplitude coefficient b A vs the interspike interval coefficient b T . In both graphs the dots represent individual parameter vectors which are coloured according to the cluster if they are fully coherent with experimental time series (selected) and grey otherwise (unselected). For each cluster in the plane (g cv , g cs ) a parameter vector has been singled out (smaller filled circles). The two larger filled circles represent the reference parameter vectors for Ca V +BK and CatSper+NHE modules. B) Illustration of the dynamics corresponding to the parameters vectors singled out in clusters c (top left), d (top right), e (bottom left) and f (bottom right) in A. The curves are the numerical solutions of the variable C in the model under the full parameter vector (coloured according to the cluster), setting g cs = 0 (black) and setting g cv = g bk = 0 (grey). Notice the overlap between coloured and black curves on the graphs on top.
show in the next section, this unique characteristic of the CatSper+NHE module cluster offers an avenue for model selection.

The CatSper+NHE module, but not the Ca V +BK module or any mixed scenario, predicts the cancelling of oscillations and sustained [Ca 2+ ] i following controlled alkalinisation
The two candidate modules for the signalling network activated by SAPs are, by construction, differently coupled to the pH i dynamics. The NHE and [H + ] i can be coupled to the Ca V +BK module in the extended model by making θ cv = 1 and θ cs = 0 in Eqs 48-52. Under these conditions the exchanger is activated by hyperpolarisation and the proton concentration becomes an output variable, in the sense that pH i does not feedback on any of the processes that determine Ca 2+ fluctuations (proton concentration H does not appear in any of the equations that govern the variables of that model variant). Therefore, this module would predict that changes in the pH i should not affect the Ca 2+ -spike trains. In contrast, in the CatSper+NHE module, the proton dynamics is tightly coupled to the Ca 2+ oscillations and is an essential part of the mechanisms underlying the oscillatory behaviour. In our analysis of the parameter dependence of the CatSper+NHE module, we systematically found oscillations in [Ca 2+ ] i that were concomitant with pH i oscillations having the same period, albeit phase differences. Parameters that lead to slower proton dynamics such that the pH i oscillations vanish, resulted in cancellation of the oscillation in [Ca 2+ ] i .
The characteristics of the two modules imply distinctive responses to a manipulation that would artificially raise and fix the pH i inside the flagellum. In this situation, the Ca V +BK module would predict that the Ca 2+ -fluctuations would be unaffected, while the CatSper+NHE module predicts that sustained pH i increase would cancel [Ca 2+ ] i fluctuations and this cation would be sustained at higher concentration. These effects are illustrated in Fig 9A and 9B, where we also overlay the results of registering Ca 2+ dynamics in the presence of ammonium chloride (NH 4 Cl), an agent that arguably raises the pH i without affecting any of the other channels modelled here. As observed in the experimental trace, the exogenous alkalinisation gives rise to a dynamics that is qualitatively similar to the one predicted by the pH i -dependent CatSper+NHE module. The conspicuous effect of alkalinisation on [Ca 2+ ] i dynamics is not predicted by any of the mixed scenarios that were otherwise coherent with the experimental time series studied the previous Sec. 2.8 ( Fig 9C). Although this result might appear surprising it has a simple rational. The oscillations in the CatSper+NHE module (as well as those obtained under all parameter vectors in cluster b) are characterised by and are strictly dependent on large amplitude alkalinisation and acidification cycles, which in turn require sufficiently long depolarisation in each cycle (Fig 9D). In contrast, the oscillations in the Ca V +BK module as well as those of the candidate mixed scenarios are characterised by very transient depolarisations resulting in sustained alkalinisation with shallow pH i oscillations. In these cases, the concomitant oscillations in the fraction of open CatSper f o cs are not driven by alkalinisationacidification cycles and therefore manipulations that keep the cytosol in an alkaline state do not block the oscillations (Fig 9D).

Coexistence of two stable states in unstimulated sperm is predicted by the CatSper+NHE module and may be revealed by transient membrane depolarisation
The CatSper+NHE module is overall the most coherent with the experimental observations. Another of its most distinctive features is that it predicts the coexistence of two possible stable equilibria characterised by [Ca 2+ ] i at either basal levels or above basal levels. This latter state is further characterised for being more acidic and with higher membrane potential than the basal state. The coexistence of two states is best seen in the bifurcation diagrams in the top graph of Fig 7. According to the CatSper+NHE module, a sufficiently strong perturbation of the membrane potential may force a switch from the lower basal steady state to the state characterised by higher [Ca 2+ ] i , where it will remain (Fig 10B). This behaviour contrasts with any model with a single resting stable state, such as the one with Ca V +BK module, in which the rise in membrane potential has, at best, transient effects on membrane potential or pH i but not on [Ca 2+ ] i , as illustrated in (Fig 10B). As a corollary of these properties, the Ca 2+ -spike train will always terminate in the basal level according to the Ca V +BK module, while according to the CatSper+NHE module, the [Ca 2+ ] i may remain at high values after the spike trains vanish

Discussion
This article explored the collection of channels that can explain the dynamics of the [Ca 2+ ] ispike trains elicited by SAP in sea urchin sperm. The large number of parameters (Table 2) used in the models studied here would deter any attempt of this kind, if it were not for the availability of quantitative time series of cGMP and membrane potential in response to several SAP concentrations. The parameters of the Upstream module that were not obtained from the literature were estimated by fitting the model to these time series data. This was possible given the structure of the signalling network in which the upstream components were assumed to be independent of the dynamics downstream and the fact that the model predicts a smooth nonoscillatory dynamics. The early signalling module was able to reproduce the long term temporal evolution of its output variable cGMP and also the early membrane potential changes in response to different SAP concentrations, measured in cell populations [10,22]. This is remarkable considering that this module makes very simple assumptions about the receptor's ligand binding and guanylate cyclase activity as well as about cGMP dynamics, neglecting adaptation mechanisms such as sensitivity readjustment of the receptor affinity that have been documented [27,28,29]. A noteworthy observation is that the cGMP dynamics at longer times, is characterised by a slow decay in which the cGMP concentration is in quasi-steady state following the slow decay of the receptor form with low guanylate cyclase activity (R L ) (Sec. 4. 2.1, Fig 2, first and second row). If the model predicted that cGMP itself would undergo periodic oscillations (arising potentially from feedbacks of pH i on phosphodiesterase activity [30]), instead of a smooth dynamics, then we would lose the inferential power of fitting the predicted cGMP value to the mean population cGMP levels.
The Upstream module alone was able to explain the early hyperpolarisation dynamics based on the joint action of KCNG and spHCN channels (Fig 2, bottom). The parameters were set such that, in the absence of any other channels, the repolarisation initiated by spHCN following the hyperpolarising current of KCNG, leads to a partial recovery of the membrane     Modelling flagellar Ca 2+ -spike trains in sea urchin sperm potential that remains below the resting value. The partial recovery implies that spHCN effective conductance density is sufficient to set into motion the downstream components (e.g. the opening of closed Ca V channels or shifting to higher opening rate of the voltage-gate of Cat-Sper) but not too strong to prevent the oscillatory [Ca 2+ ] i dynamics. This property of the signalling network is amenable to experimental corroboration by assessing whether the membrane potential recovery is only partial under a setting in which all calcium entry is prevented by using zero external concentrations of this cation. The upstream SAP-induced cGMP response curve is linked to the downstream [Ca 2+ ] i fluctuations by the dynamics of the network. The key question was: what components of the downstream modules correctly predict this link? Our analysis indicates that either one of two modules featuring voltage-dependent channels could be responsible for the dynamics of [Ca 2+ ] i fluctuations. The first module calls into action Ca V channels that undergo voltagedependent sequential transitions from inactivated to closed, to open and back to inactivated states. These irreversible transitions bring forth a limit cycle with constant period with an amplitude that is controlled by cGMP activity. These Ca V channels, which are predominantly inactivated at the resting potential, are released from inactivation upon a transient hyperpolarisation and subsequently opened when spHCN channels activity raises the membrane potential. In this module, the oscillations are strongly modulated in amplitude and tempo by the Ca 2+ -dependent channel BK to produce the experimentally observed progressive increase in the intervals between [Ca 2+ ] i spikes. The oscillations predicted by this module vanish when cGMP falls below a critical value involving a Hopf-type bifurcation (Fig 5A).
The second module capable of explaining the organisation of the [Ca 2+ ] i fluctuations features a single channel, CatSper, that is voltage-and pH-dependent and is reversibly inactivated by intracellular Ca 2+ ions. In contrast with the previous module, which is set in motion directly by the effect of transient hyperpolarisation that relieves Ca V channels from inactivation, the initial opening of CatSper channels is triggered indirectly by alkalinisation of the cytoplasm driven by the hyperpolarisation-dependent Na + /H + -exchanger activity (S3 Fig). The cycling behaviour is produced by the negative feedback that depolarisation driven by CatSper activity has on inactivating the Na + /H + -exchanger, reducing pH i and decreasing CatSper activity. Another feedback arising from the inactivation of CatSper by [Ca 2+ ] i plays a more preponderant role in controlling the length of the Ca 2+ -spike train and in increasing the interspike intervals. As in the case of the first module, the oscillations generated by the module stop when the cGMP activity falls below a critical value (Fig 7).
It is noteworthy that, under the reference parameters used in this article, the overall temporal scale of the calcium response produced by the CatSper+NHE module is slower compared to the experimental tracings, in terms of real time units; however, if the relative progressive increase of interspike intervals is contrasted instead, this property becomes more marked in this model than in the one featuring the Ca V + BK module, and more closely related to the experimentally observed (Fig 3C). In the numerical solutions with the former module the interval between the last spikes can be two or three fold larger than the interval between the first spikes in the train, as observed in the experiments in single cells. A comparable increase in interspike intervals can be obtained with the Ca V +BK module under different parameter vectors (such as those in cluster a). Conversely, the progressive decrease in the amplitude observed experimentally is better accounted for by the Ca V +BK module than by the CatSper+NHE module under all the parameters studied here.
These two downstream modules confer some distinguishable properties to the signalling network when coupled to the Upstream module. One of those properties is the envelope of the membrane potential oscillations. In the module featuring the Ca V + BK channels pair. This envelope has a shape similar to that of the [Ca 2+ ] i spikes. In contrast, the CatSper+NHE module predicts that the amplitude in the membrane potential peaks should increase until it suddenly ceases (Fig 6, bottom), thus establishing an opposite trend to the calcium-spike train envelope. Interestingly, the ascending V envelope is more similar to the observed at prolonged times in membrane potential measurements of sperm populations stimulated with high doses of SAP [10].
Another good example of distinctive features between the alternative scenarios is the permissive cGMP range that allows for oscillatory solutions. At too low (close to resting) or too high intracellular GMP concentrations, the model with CatSper predicts a sustained rise [Ca 2+ ] i without noticeable oscillations (Fig 7). Conversely, in the Ca V + BK module, the periodic regime is maintained for any value above a critical level of cGMP, however, above more elevated concentrations, the amplitude and period saturate to nearly constant values (Fig 5A). These properties of the model could be addressed experimentally by phosphodiesterase inhibitors or by elevating cGMP to saturating concentrations by uncaging an analogue. It is worth recalling the reports that intracellular uncaging of cGMP leads to an increase in [Ca 2+ ] i that is oscillatory in the case of A. punctulata [9,31] but transient in the case of S. purpuratus ( [15]) spermatozoa. The discrepancy between the two species could be reconciled if there are speciesspecific differences in the boundaries of cGMP intervals that produce oscillations or alternatively by experimental differences in the uncaged concentration of the cyclic nucleotide.
The difference between the two alternative network modules that is most decisive is the distinct interplay with cytosolic pH and proton dynamics.While the periodic alkalinisation and acidification is an essential part of the limit cycle characteristic of the CatSper+NHE module network, proton dynamics is secondary and inconsequential in all the networks that included the alternative Ca V + BK module. For this reason, only the CatSper+NHE module predicted that externally forced cytosolic alkalinisation should lead to sustained [Ca 2+ ] i as observed experimentally following the addition of NH 4 Cl to S. purpuratus spermatozoa (Fig 9A). This decisive results prompts us to conclude that the channel responsible for the [Ca 2+ ] i responses to SAP is the pH-dependent CatSper instead of Ca V . However, this conclusion, does not follow simply from the pH sensitivity of CatSper or from the SAP induced alkalinisation. These features were also present in the mixed scenarios that failed to reproduce the response to forced alkalinisation, despite being compatible with the other observations. Further considerations on this conclusion are called for. One of its corollaries is that observing [Ca 2+ ] i oscillations dependent on CatSper and NHE implies necessarily the concomitant observations on pH i oscillations with the same period. Oscillations of pH i with the same period of those of [Ca 2+ ] i were not detected in experiments where the levels of pH i have been monitored in individual cells [32]. This might reflect limitations of the experimental system (e.g. too low signal to noise ratio of the pH probe as compared to the Ca 2+ probe, particularly if the relevant proton concentration is the one near the membrane and not in the bulk of the flagellum) or, alternatively, it might mean that intraflagellar pH is not oscillating periodically. In the latter case, the hypothesis that CatSper is the main responsible for fluctuating Ca 2+ currents, at least in the way we modelled it here, is questionable and must be reexamined opening a new avenue for modelling research. In this context, it might be relevant to note that oscillatory fluxes of Ca 2+ and protons have been reported with the same period as growth oscillations of pollen tubes (reviewed in [33]).
The observation that a rise in pH i precedes the first Ca 2+ spike in sea urchin sperm stimulated by SAP has been interpreted as an indication that Ca 2+ currents are a consequence of this early alkalinisation [34,12,32]. This interpretation became a cornerstone of the hypothesis that pH-sensitive CatSper channel is responsible for Ca 2+ fluctuations. It is worth calling attention to the early time courses of pH i and [Ca 2+ ] i obtained by the numerical solution of the model featuring the Ca V +BK module (Fig 9A). In this model, where Ca 2+ currents are, by construction, independent and unaffected by pH i one can nevertheless observe an onset of alkalinisation that precedes the first [Ca 2+ ] i spike. This illustrates the frailty of inferring causation from time series data in general and emphasises the need to support the hypothesis that a pH-sensitive Ca 2+ channel is involved responses to SAP with more decisive experiments such as the controlled pH i manipulation, as proposed here.
Another key prediction of the CatSper+NHE module is the existence of bistability at the resting values of cGMP. This implies that perturbations exist that may lead to a sustained [Ca 2+ ] i value above the normal resting one, in the complete absence of SAP signals. We illustrated this by a transient rise in membrane potential that locked the system in the stable state characterised by higher, tonic-like [Ca 2+ ] i (Fig 10B top). Another perturbation with comparable effects would be a transient yet complete depletion of cGMP that would push the system across a bifurcation point in which the lower equilibrium vanishes. Following such perturbation, hysteresis is predicted such that the system could remain in the higher [Ca 2+ ] i state after restoration of the basal cGMP levels. The CatSper+NHE module bistability also implies that after a [Ca 2+ ] i -spike train terminates, the system may or may not return to the basal level. It is interesting to note that this alternative [Ca 2+ ] i equilibrium state also entails low pH i levels and depolarised membrane potential (Fig 10B bottom). If such a state existed, its physiological meaning would be associated with a state of quiescence, since the acidification of the cytoplasm below the basal pH i inhibits both the metabolism and the dyneins, which are the molecular motors that drive the flagellar beating [35,36]. On the other hand, the viability of sperm would be affected in this state, since sustained high levels of calcium represent a universal signal that triggers cell death in the majority of eukaryotic cells [37,38,39]. These properties of bistability and hysteresis at cGMP levels close to the resting state are not present in the Ca V +BK module, thus offering yet another way of teasing apart the two hypotheses.
Our analysis indicates that proportionate mixed scenarios involving components of the two modules are implausible in the sea urchin sperm. The Ca V and BK channels can fine-tune the temporal structure of spike trains elicited by CatSper if present at minute densities. At densities or currents comparable to those of CatSper these channels will disrupt the [Ca 2+ ] i oscillation train or bring forth an oscillatory regime in which the proton dynamics plays an irrelevant subsidiary role. In the same vein, the CatSper+NHE module predicts [Ca 2+ ] i -spike trains in agreement with the experimental observations only if BK channel densities are negligible. If one assumes that the unitary conductance of BK channels in marine invertebrates is closer to the estimated in mollusk and crustacean cells [40,41,42] (under asymmetric conditions that approach physiological gradients) than its counterpart in mammalian cells [43], this would be of the order of �70 pS. Under these conditions, the BK conductance densities that allow for oscillatory trains would imply channel counts of less than 1 molecule per flagellum, given the membrane area (Table 2). Therefore, although BK would play a key role in a scenario wherein Ca V channels were the ones driving calcium, the role of this channel, if any, should be minimal if CatSper turns out to be the predominant channel. The role for BK in [Ca 2+ ] i responses to SAP was inferred indirectly by the alterations of the period of these oscillations by blockers such as niflumic acid [11,44] and Iberiotoxin [17], which were corroborated by discrete logical network models. In addition, proteomic studies in sea urchin sperm have detected BK [18]. In the present quantitative framework, there are shared effects of BK both on the Ca V +BK and CatSper+NHE modules, such as defining a tonic-like minimal value for the nadir of [Ca 2+ ] i oscillations and increasing the period of the oscillations. This notwithstanding, the CatSper+NHE module precludes the presence of BK in sea urchin sperm at significant levels.
The use of models based on differential equations allowed us to confront their predictions with several quantitative data sets, gaining insights beyond those stemming from logical network models. The slow quantitative dynamics in cGMP as well as the trends in the spike amplitudes and interspike intervals of the [Ca 2+ ] i -spike trains, which were instrumental in our analysis, are a kind of rich information source to which the discrete logical network modelling is blind. We also improved the continuous model of the [Ca 2+ ] i core oscillator in urchin sperm developed before [21], not only by bringing the network composition to the state of the art knowledge, but by grounding the analysis on experimental measurements at multiple levels of the network. Despite these new insights provided by ordinary differential equations, we foresee future directions for theoretical studies using other quantitative frameworks. The reports on the spatial localisation of CatSper channels in 4 distinct bands along the flagellum of both mouse [45] and human [46] sperm is fascinating and calls for partial differential equations to understand the significance of such topological organisation. In our models we assumed that molecular species such as cGMP, Ca 2+ ions and protons were diluted uniformly within the volume of the flagellum. This likely imposed quantitative constraints on parameters that might be alleviated by considering local concentration asymmetries arising from the spatial arrangement of the membrane channels relative to each other (e.g. in human sperm Hv1 channels, the putative major pH i regulator, are located between two of the CatSper bands [46]). In restricted domains, ample oscillations in ion concentrations may be enabled by smaller conductance densities than the ones required by our mean field model. Whether such a spatial arrangement of multiple channels exists in sea urchin sperm remains an intriguing question. The variance in the timing and amplitude of [Ca 2+ ] i spikes in individual cells, which was neglected in our analysis of the overall trends of the spike trains, may contain information on the underlying stochastic processes that could be better captured by Gillespie-type stochastic simulations or spatial-stochastic models.
We used bifurcation analysis of a subset of the differential equations of the full system as a heuristic tool to understand how the temporal evolution of cGMP activity affects the envelope and intervals between [Ca 2+ ] i spikes. It is important to bear in mind that the present analysis involves an asymptotic approximation, i.e. in our case G was set to a fixed value as a way to transform this variable into a parameter and quantify the amplitude and period of the stable limit cycles. The experimental corroboration of such analysis would require development of a cGMP-clamping technique, analogously to classic voltage-clamping experiments or synthetic biology approaches that seem presently unfeasible in sea urchin sperm.
We end by noting that the present work illustrates the added value of quantitative analysis which beyond offering a framework to characterise mechanisms overcoming experimental limitations, helps to identify the main elements controlling the system. Quantitative modelling allows to define guidelines for future experiments to assess the theoretical predictions, to identify key properties and testable quantities, and to prioritise these properties and quantities according to their impact and significance for the system.

The general modelling framework
Our strategy is to describe the dynamics of the SAP-activated signalling pathway in order to understand how the [Ca 2+ ] i fluctuations are controlled by different molecular components. A system of ordinary differential equations (ODE) describes the dynamics of the membrane potential and of the intraflagellar concentrations of cGMP, protons and Ca 2+ . These four essential variables depend on the activities of different flagellar channels and enzymes, the dynamics of which are also described by ODE (Sec. 4.2). The four variables are the main observables in the model and offer the means to assess its predictions by comparison with experimental data.

Membrane potential.
The dynamics of the membrane potential, V, is described with a Hodgkin-Huxley formalism [47]. The temporal derivative of V is defined as a sum of ion current densities normalised by the flagellum specific capacitance, according to Kirchoff's law of charge conservation: Each current density, indexed by i, is associated to a given channel type to be chosen from a set of ion channels, i.e. i 2 {kn, kc, cv, bk, cs}. The current terms I i are defined according to Ohm's law: where the E i is the reversal potential of the channel i, g i is the channel-specific maximal conductance density defined as the product of the channel unitary conductance by the channel density in the flagellum membrane, and f o i is the fraction of open channels conducting current. The aggregate parameter g i is a convenient description of the conductance density given that for most channels in sea urchin spermatozoa, the unitary channel conductance, the channel density or both are unknown. The fraction f o i is a variable whose dynamics reflects the complex regulation of channel gating. Eq 2 is therefore a general framework and channel specific details are presented individually in Sec. 4.2.
Following the Hodgkin-Huxley tradition, we define a leakage current term I L = g L (V − E L ) that ensures that the resting membrane potentialV ¼ E m is asymptotically reached by imposing the following constraint on E L , given the parameters and variables f i at equilibrium: The leakage current is interpreted as the net current arising from all the channels and transporters that are not explicitly described in the model. As a remark on notation, in this article a hat-decorated variablex refers to the equilibrium value of the variable x.

Cyclic nucleotide concentration.
The concentrations of the nucleotides cGMP and cAMP increase in the sperm flagellum in response to SAP stimuli. These nucleotides affect the permeability of several nucleotide-gated ion channels, in turn changing membrane potential. Considering that cGMP is the nucleotide showing the highest changes after SAP binding and that the experimental release of its analogues in the flagellum is sufficient to produce a dose response that closely mimics the response elicited by the physiologic ligand, we restricted our analysis to this nucleotide, neglecting the eventual role of cAMP.
The intraflagellar concentration of cGMP is denoted by G and its dynamics described by the following differential equation: In this equation, it is assumed that the activities of general guanylate cyclases and phosphodiasterases lead to a constitutive turnover of cGMP, with a constant source σ G and turnover rate δ G . In the absence of SAP signals, this basal turnover results in a cGMP at rest stationary concentration G r = σ G /δ G . The basal turnover is assumed to be perturbed by the additional production of cGMP due to active SAP receptors, described by the term ρ G (Sec. 4.2.1). A phosphodiesterase that impinges on sperm motility has been described in sea urchin [30]. Although its GAF domains could indicate allosteric regulation of catalytic activity by cGMP itself, the phosphodiesterase response to cGMP gave no hint of such nonlinearity, justifying the assumption of constant turnover rate δ G . The same report has shown that the potential regulation of the phosphodiesterase activity by pH [30], and this activity could be also regulated by [Ca 2+ ] i via calmodolin, as documented in other cellular systems and species [23]. Eq 4 neglects any potential effects of pH i and [Ca 2+ ] i on phosphodiesterases, which is a simplifying assumption that proved to be reasonable by comparison with quantitative data.
The variable G was compared with measurements of total cGMP by radioimmunoassay in bulk cell experiments previously reported [22]. Since G represents the effective intraflagellar concentration of cGMP available to bind to KCNG channels, it was necessary to introduce a scaling factor for the term ρ G . This factor, θ G , consists of the product of flagellar volume (v f ), cGMP buffering capacity in the flagellum (B G ) and Avogadro constant (N A ). In this way, only a fraction of the total cGMP synthesised by the active receptors participates in the signaling response and the G signal is kept within the nM sensitivity range reported for KCNG channels [48] (Sec. 4.2.2).

Proton concentration.
The permeability of some channels is affected by cytoplasmic concentration of protons. The intraflagellar proton concentration, denoted by H, is described by the following differential equation: which assumes that there is a source of intraflagellar protons, σ H , a basal efflux rate δ H and a variable efflux J H mediated by the activity of the sodium/proton exchanger (Sec. 4.2.4). The parameter values were constrained by imposing the basal proton concentration to be 7.9 × 10 -8 M (pH i = 7.1, Table 2), obtained by solving H for equilibrium. The relative changes of the variable H in time were further compared with pH-sensitive, fluorescent probe signals measured in bulk populations [34] or in individual cells [32].

Intraflagellar Ca 2+ concentration.
The intraflagellar Ca 2+ concentration, represented by the variable C, is the main observable of our study. Its dynamics is described by the following differential equation: with j 2 {cv, cs}. Just as in the case of the previous components, we assume that different sources and sinks of intraflagellar Ca 2+ result in a net constant influx rate σ C and a basal efflux rate δ C . The [Ca 2+ ] i dynamics is further controlled by the activities of specific Ca 2+ channels, described by the third term in the equation as the sum of ionic current densities scaled by the ratio of flagellum membrane area over the product of ion charge of a Ca 2+ ion, flagellar volume and Faraday constant. The list of currents to be considered will depend on the scenario being analysed.
The variable C can be compared to measurements by Ca 2+ -sensitive fluorescent probes, either in individual cells or in populations [14,9,22]. Because the population measurements give only information about the trend of the envelope of the spike train, we relied on individual cell measurements for the analysis of the temporal structure and the amplitude of the Ca 2+ spikes, as described in Sec. 4.4.

Signalling components and channels
In the following section, we will present the molecular components and channels that were implicated in SAP signal transduction and that are relevant for our study. For each of these components we will describe: i) the evidence for their role in SAP-signalling, ii) a mathematical model for its dynamics, iii) the reference parameter values and how they were obtained. 4.2.1 SAP-receptor and guanylate cyclase. The binding kinetics of SAP to its respective sperm receptor has been studied in sea urchin species such as L. pictus [49] and S. purpuratus [50], whose spermatozoa respond to the decapeptide speract. The reported values for the association and dissociation rate constants for speract imply that the formation of a ligand-receptor complex is essentially irreversible during a chemotactic response, with an average life time of several hours. Since the time scale of the SAP response is of the order of seconds, for the sake of simplicity we neglected the dissociation rate. In order to explain the kinetic features of experimental cGMP data [22] we propose the following mechanism: once free SAP molecules S bind to free receptors R F , they go through several states in an irreversible manner, each with decreasing guanylate cyclase activity, namely high activity R H , low activity R L and inactive R I receptor forms. In the context of previous work, the molecular species R H and R L might be considered as equivalent to the phosphorylated and dephosphorylated forms of guanylate cyclase, respectively [51,52,53,28]. Receptor activation mechanism is represented by the following kinetic diagram: which is translated in the following ordinary differential equations: subject to the mass conservation relations: with S being the concentration of free SAP S and R F , R H and R L and R I being the number of the R F , R H , R L and R I receptor forms per cell, respectively. R T and S T are constants denoting the total receptors per cell and SAP concentration, respectively. As in other studies of soluble ligand-membrane receptor interactions [54], we describe the amounts of soluble SAP molecules S in molar concentration and those of the membrane receptor forms, R T , R F , R H , R L and R I , in molecules per cell. As the SAP binds to the membrane receptors, its concentration in solution decreases proportionally to the bound receptor forms with a constant θ R , and this decrease is described by the second term of the conservation Eq 13, θ R (R H + R L + R I ). This representation is convenient to model bulk studies as we set θ R = s/N A , where s is the number of sperms per litre and N A is the Avogadro constant.
Regarding the activation of guanylate cyclases stimulated by SAP, minor differences exist at the molecular level between sea urchin taxonomic orders. In sperm of S. purpuratus and L. pictus, which belong to the Echinoida order, receptor and guanylate cyclase have been reported to be two separate yet presumably tightly coupled membrane-bound proteins; once a SAP molecule binds to the former [55] this in turn stimulates the activity of the latter [56]. In the third species, which belongs to Arbacioida order, a single protein has the activities of ligand binding and guanylate cyclase [57,58]. To capture the commonalities in a single model, we only take the latter and simpler case, thus ignoring any possible delay between SAP binding and guanylate cyclase activation. Thus, the term corresponding to the SAP-dependent cGMP synthesis, σ G , in Eq 4 is defined as: where the constants k H and k L measure the guanylate cyclase activities of the high (R H ) and low activity (R L ) SAP-receptor forms, respectively.

Cyclic nucleotide gated K + channel.
The K + -selective cyclic nucleotide gated (KCNG) channel is a pseudotetrameric channel with up to 4 cyclic nucleotide binding domains [59,48], nonetheless, the binding of only a single cGMP molecule is necessary and sufficient to open a channel, thus giving a non-cooperative gating unlike most cyclic nucleotide gated channels [60,48]. We therefore assume that the fraction of open channels, denoted as f o kn , increases with a rate directly proportional to the cGMP activity (G) and decreases with first order kinetics.
Expression for both the steady state fraction of open channels,f o kn , and the characteristic time, τ kn , can be derived from the opening and closing rates: with K kn ¼ b kn a kn . The ion current density through these channels is modelled as: where E K is the Nernst potential of the K + and g kn denotes the maximal conductance density of the KCNG channel (the product of the unitary conductance by the density of channels).

Hyperpolarisation-activated cyclic nucleotide-gated channels.
Hyperpolarisationactivated cyclic nucleotide-gated channels from sea urchin sperm (spHCN) are cationic channels, weakly selective and carry a mixed current of Na + and K + with a permeability ratio P K / P Na of *4-5 and a reversal potential of E hc = -30 mV-−10 mV which would allow an inward Na + -current under physiological conditions that, when open, depolarises the cell membrane [61,62,63]. Unlike typical spHCN channels, spHCN-current inactivates in the absence of cAMP, and the latter increases the maximal open-probability instead of shifting the half-activation voltage of the G/V curve [62,64]. Notwithstanding that, for the ionic current equation associated with this channel, we assumed in the model that cAMP is always present at sufficient concentration to avoid current inactivation, hence any role of cAMP in the overall dynamics is neglected. Following a classical Hodgkin-Huxley type model, we propose that spHCN opening depends on three independent gating variables, denoted by m hc , such that the fraction of open channels is f o hc ¼ m 3 hc . The dynamics of the gating variables is described by: wherein the steady state fraction of active gating variable,m hc , is modelled with a Boltzmann equation, whereas the characteristic time, τ hc , is modelled as a Gaussian function of the voltage:m The current via spHCN channels I hc is given by: 4.24 Electroneutral sodium/proton exchanger. A sperm-specific electroneutral sodium/ proton exchanger (NHE) is expressed in sperm flagella membrane [65]. Its activity is enhanced with hyperpolarised membrane potentials and inhibited with depolarisation. When active, it makes use of the Na + concentration gradient between external and internal media to extrude protons by means of a Na + /H + electroneutral exchange mechanism with a 1:1 stoichiometry, consequently increasing pH i [66]. Even though NHE is not a channel, it is distinguished from the rest of NHE family members in having a voltage sensor domain that is homologous to that of voltage-gated ion channels [13], which might explain the observed voltage-dependent Na + /H + exchange in sea urchin sperm in spite of being an electroneutral net exchange [67,66,68,69]. We therefore model this voltage dependency as a regular gating mechanism similar to that of Hodgkin-Huxley like equations. The fraction of active NHE molecules, denoted f a nh , is described by: where the activation and inactivation rates are the following functions of the membrane potential: Proton flux as a consequence of exchange activity of NHE, designated J H (mol s −1 ), is modelled with reversible, rapid-equilibrium, random order bi-bi kinetics [70]: where: J max (mol s −1 ) is the maximal exchange rate; [Na + ] o and [Na + ] i are the extracellular and the intraflagellar concentrations of Na + , respectively; pH o is the extracellular pH; and K Na is the dissociation equilibrium constant for Na + .

Voltage-gated calcium channels.
There has been suggestive evidence for the involvement of classical voltage-gated calcium channels in the SAP-activated response, i.e. Low-Voltage-Activated [14,10] and High-Voltage-Activated channels [16]. Here we focus on T-type calcium channels, which belong to the first class and are known by being associated to spiking activity in other cell types. We model their kinetics by the following diagram, in which illustrates that they transit irreversibly through three states, namely, inactive Ca The gating rates are themselves voltage-dependent exponential functions similar to those of Hodgkin-Huxley (HH) type models: Denoting the fraction of open, closed and inactive channels as f o cv , f c cv and f i cv , respectively, we have two differential equations for the kinetics of the first two states, whereas the last one is calculated from the conservation equation After solving gating equations for equilibrium, we get: The ion current density through these channels is modelled as: 4.2.6 CatSper channels. The sea urchin genome contains the gene set necessary for the expression of the four core and the three auxiliary subunits that constitute the CatSper channel [71]. Furthermore, the expression of these subunits has been demonstrated in A. punctulata sperm flagella [12] and most of them in S. purpuratus [18]. In order to model CatSper gating, we consider three regulatory mechanisms reported in studies from mouse, human and sea urchin sperm: voltage-dependence, modulation of voltage-sensitivity by pH i , and inactivation by intracellular Ca 2+ [72,73,74,12]. We implement these mechanisms using two independent gating variables, m cs and h cs , representing the fraction of channels with their gate modulated by both voltage and pH i open, and the fraction of channels with Ca 2+ -dependent gate that is not inactivated by calcium. The fraction of open CatSper channels conducting Ca 2+ currents, f o cs , is given by the product of these gating variables: whose dynamics is described by the following differential equations: The opening and closing rates of the voltage-and pH-dependent gate were chosen to be the following symmetric functions: such that the equilibrium open gate curve is a Boltzmann function of the voltage: The dependence on pH i is ensured by defining the half-maximal voltage, v cs , as the following function of the proton concentration H: v cs ¼ s 3 s s 5 4 ðÀ log 10 HÞ s 5 þ s in order to produce a shift of the gating curve towards more negative membrane potentials upon alkalinisation of cytosol, as observed in [72,73,12]. The second gate, inhibited by [Ca 2+ ] i , is assumed to switch between open and closes states with the following rates: b h ¼ s 10 ðC=s 7 Þ s 9 ð42Þ whose function forms were chosen such that when solving for the open gate fraction at equilibrium,ĥ cs , one obtains a Hill-function: The ion current density through CatSper is modelled as: 4.2.7 Big-conductance potassium channels. Big-conductance potassium channels (BK, also known as Slo1) are activated by [Ca 2+ ] i . These channels have been implicated in sea urchin sperm signalling by pharmacological blocking assays [11,17] and adopted in discrete models [17,20,18]. For their dynamics we assume two states, open or closed, where the open-channel fraction, denoted by f o bk , is described by the following differential equation: with the closing rate kept as a constant β bk and the opening rate being the following Ca 2 + -dependent function: The ion current density through these channels is expressed as:

Computational implementation and parameter estimations
The systems of differential equations were solved numerically with Mathematica 11 (Wolfram Research, Inc., Mathematica, Version 11.0, Champaign, IL), XPPAUT 8 [24] or R (using the package deSolve [75]). XPPAUT 8 was used for the bifurcation analysis (see below). Mathematica was further used to fit the model to experimental data using the Nonlinearmodelfit function with the optimisation method Differential Evolution with constraints. In the analysis, we defined several model variants based on subsets of the variables and equations presented above. Some of these defined modules were used to estimate the parameters and some were used to analyse specific properties. The Upstream module (illustrated in Fig 1) was fitted to data in two phases: • The state space comprised by the variables {S, R F , R H , R L , G}, which includes the reactions involving the SAP-receptor (Eqs 8-11) and the cGMP dynamics (Eq 4), was fitted to kinetic data of cGMP measured in bulk cell populations [22]. Data points were read from the indicated publication and are expressed in units of pmol per 10 8 cells. The cGMP kinetics obtained with the different SAP concentrations used in [22] were fitted simultaneously with the model (S1 Fig). During this fitting the dimension cGMP was pmol per 10 8 cells instead of molar concentration, such that the scaling factor θ G (described in the Sec. 4.1.2) was not called in. The parameter value of σ G obtained in this phase, being measured in pmol per 10 8 cells per s, was subsequently rescaled to molar concentrations per s, taking into account the flagellar volume.
• The remaining variables of the module (f o kn , f o hc , V), which involve the early electrophysiological changes caused by SAP (i.e. hyperpolarisation and repolarisation due to the opening of KCNG and spHCN, respectively), were used to fit the gating parameters of the ion channels. The spHCN parameters were either extracted or calculated from [62,63,61,76], while those of KCNG were based on the values reported in [10,59,48]. Fixing these values, the conductance density parameters of these channels were calibrated together with B G (included as a free parameter in the scaling factor θ G ), by fitting the dynamics of V with the kinetic data from [10], which comprise the early response of the membrane potential change. Only the first second of the time series was used in the fitting, taking into consideration that, by itself, this module is expected to fit only the early transient behaviour of V before other downstream channels with strong effect come into play.
In order to assess and calibrate this module, we assumed that the overall dynamics under both single-cell and population regimes can be compared using a scaling factor; this depends on the flagellar volume term to convert the cGMP signal into the effective concentration capable of exerting a physiological effect at the single cell level (see parameter θ G , Sec. 4.1.2).
To study the dynamics of intraflagellar calcium, we set out to explore the possible contribution of different channels under two main scenarios.
• In one scenario, the model features the Ca V +BK module and has the state space {S, The values of the Ca V channel parameters were adjusted based on the characteristics of low-voltage-activated T-type Ca 2+ channels [77], for which there is evidence suggesting their presence in sperm [14,10] (Table 1)

Comparison of the structure of the Ca-spike trains in experiments and models
The [Ca 2+ ] i dynamics in response to SAP was imaged in individual S. purpuratus sperm immobilised on a cover slip using intracellular Ca 2+ -sensitive fluorescent probe following the uncaging of a photoactivatable Speract analogue by a 0.3s pulse of UV light (experimental details in [78]). By measuring the average fluorescence of the flagellum of each cell at each time point we obtained a time series of this proxy of [Ca 2+ ] i ). Each time series was smoothed using a cubic spline (as implemented by the R function smooth.spline) and the maxima of the smooth curve used to identify Ca-spike trains that were characterised as (t k , A k , T k ), where k is the ordered spike index (k = 1, 2, 3, . . .) and t k , A k , and T k = t k+1 − t k are the time, the amplitude and interspike interval. To summarise the structure of each individual time series we computed the average interspike time T k , and the regression coefficients of the relative amplitude A k /A 1 and relative interspike interval T k /T 1 over the spike index, denoted b A and b T , respectively. To compare the [Ca 2+ ] i dynamics in the model the same procedure was applied to the numerical solutions of C. The structure of the spike trains in the model was considered to be coherent with those experimentally observed according to the following selection criteria: the average spike interval is in the range 0.37s-1.38 s, the predicted (b T , b A ) point falls within the convex hull defined by experimental values in the bottom right quadrant, i.e. the spike amplitude decrease whereas the interspike interval increase in time b A < 0 and b T > 0. As mentioned in the Sec. 2.3, this pattern of the spike trains is the most representative of the experimental time series and the most challenging for the models. In addition, to filter out spontaneous oscillations and alternative equilibria in multistable regimes, we imposed that the resting steady state value of [Ca 2+ ] i ,Ĉ, obtained following numerical solutions with S = 0, is about 100nM. It is important to notice that the fluorescent probe intensities were not calibrated and therefore mean flagellar intensities are not accurate [Ca 2+ ] i measurements. We assumed that the fluorescence intensity was proportional to the actual concentration of [Ca 2+ ] i , focusing our analyses on relative differences and values. This assumption underlies the co-plot time series and model solutions presented in the figures and also the direct comparison of the theoretical and experimental b A values. We believe this assumption is reasonable, considering that all Ca-fluorescence intensity spikes analysed were below the control intensities elicited by ionomycin.

Mixed scenarios parametrisation and parameter sensitivity
We analysed mixtures of Ca V +BK and CatSper+NHE modules by assembling the most extended model with state space {S, The parameters of this model are the union of all the parameters listed in Table 2.
Any mixed scenario involving the two modules can be parameterised by setting: g bk ¼ y cv g cv bk ð49Þ g cs ¼ y cs g cs cs ð50Þ g L ¼ y cv g cv L þ y cs g cs with the superscripts cv and cs denoting the reference values in the Ca V +BK and the CatSper +NHE modules, respectively, and all the remaining parameters taking the values in Table 2. θ cv and θ cs are non-dimensional scaling factors that allow to weigh the contribution of each module. With this formulation, the original Ca V +BK module dynamics can be recovered in this extended model by setting θ cv = 1 and θ cs = 0, whereas the CatSper+NHE module can be obtained with θ cv = 0 and θ cs = 1. From the above equations, we can express a ratio of the maximal conductance density of CatSper with respect to that of Ca V : In a first approach to mixed scenarios we interpolated between the two modules. Starting with the CatSper+NHE module we added fractions of Ca V and BK channels in three ways: a) adding Ca V + BK while decreasing proportionally CatSper by setting θ cv = θ and θ cs = (1 − θ) with θ 2 [0, 1] (S4 Fig); b) adding only Ca V by setting θ cv = θ and θ cs = (1 − θ) while forcing g bk = 0 (S5A Fig); and c) adding only BK by defining g bk ¼ yg cv bk and setting θ cv = 0 and θ cs = 1 (S5B Fig). A more systematic analysis of mixed scenarios contemplated the conductance of channels that compose the two modules but also variants of the values of all the other parameters, except for those representing the Upstream Module (which were fixed at their point estimates in Table 2).The exploration of this parameter hyperspace was done using an algorithm inspired on adaptive evolutionary dynamics, specified in S1 Algorithm. Briefly, the algorithm comprehends two phases, one of neutral exploration of the parameter space by evolving a seed population of parameter vectors in a regular grid, followed by further exploitation of parameter vectors selected according to the criteria described above (Sec. 4.4). First, each parameter vectorl in the nodes of the grid was obtained by taking θ cv , θ cs 2 {0.0001, 0.1, 0.2, 0.3, . . ., 2.0}. Each vector in this population was used randomly as a seed to generate a daughter parameter vector obtained by multiplying each parameter value by an independent lognormal distributed random variable (specified as e x with x being normally distributed with mean 0 and standard deviation 0.04). Each daughter parameter vector was appended to the evolving population, irrespective of the output behaviour, serving as random seed to generate new daughters. The method to generate parameter vector variants favoured values close to the reference ones (Table 2), reducing the probability that the random exploration would spread to unrealistic regions of the parameter space (e.g. leading to inactivation, opening and closing of Ca V within incongruous membrane potential intervals). The dynamics predicted by each parameter vector l, seeds and daughters alike, was obtained by solving numerically the extended model. The process of generation of new parameter vectors was repeated until the clouds of points, spreading from the nodes of the grid, once projected in the plane (g cv , g cs ), started to visibly coalesce. At this stage, several parameter vectors satisfied the selection criteria for coherence with experimental observations and were selected for further exploration in the second phase of the algorithm. In each iteration, a vector in the selected population was randomly picked to generate a daughter vector by the same procedure as above. This daughter vector was analysed and added to the evolving selected population, serving as seed for further exploration, if it met the selection criteria; otherwise it was stored in a list of unselected parameter vectors. The latter vectors are by construction closely located in the parameter space to those of the evolving selected population from which they derived but, being incoherent with observations, they were not further explored. The selected and unselected vector populations were used to build clusters based on the proximity in the projection plane (g cv , g cs ). The iterative generation of daughter parameters from the selected population was repeated until the clustering in the plane (g cv , g cs ) was stable (i.e. no new clusters were observed after a large number of iterations). Within each final cluster, the fraction of number of selected parameter vectors over the total number of vectors is a kind of measure of the parameter sensitivity. The smaller this fraction the more sensitive the model is within the cluster.

Resting state and initial conditions
Initial conditions of the gating variables of KCNG, spHCN, Ca V , BK and NHE, as well as parameters σ C , E L and δ H , were all calculated by solving the system at equilibrium in the absence of SAP, that is, setting each differential equation to 0 under the condition S = 0. Physiological observables G, C, H and V were initialised to their respective resting values, as reported in the literature ( Table 2).
The parameters governing NHE activation were chosen such that approximately 40% of exchangers were active at the resting state, as inferred from experimental data in activated sperm without any chemoattractant stimulation [35,79], whereas the maximum effective proton flux by NHE, J max , was fitted to limit the pH i up to *7.7, after physiological stimulation by SAP [80,81].  [22] and the lines are the single cell numerical solutions of the variable G in the Upstream module, with state space {S, R F , R H , R L , G} that best fit the ensemble of the data (i.e. the four curves simultaneously). The values of R T , r 2 , r 3 , δ G , k L and k H were fitted while fixing the remaining parameters (Table 2). Because the measurements of the SAP-receptor association rate constant r 1 are technically more reliable than those of the number of receptors R T (reported values include 14000 [88,22], 3 × 10 5 [28] and 1 × 10 6 [27] molecules per cell), we fixed the value of r 1 and fitted R T . We fixed also the value of the factor that converts the amount receptor-bound SAP molecules to soluble concentration θ R = s/N A = 6.14 × 10 −5 nMcell based on the reported cell density s = 3.7 × 10 7 cellmL −1 [22]. Finally, G r was preset to 1.24 nM [22] and the condition σ G = G r /δ G held during the fitting. (TIF) S2 Fig. Calibration of spHCN gating parameters. In the upper panel, experimental traces of ionic currents measured by whole cell patch clamp technique in HEK cells expressing heterologously spHCN and loaded with photoactivatable cAMP analog, which in turn was uncaged by UV light. The set of currents correspond to different voltage pulses (indicated at the end of the trace). Data extracted from figure 4a of [62] (gray lines). Each trace was fitted to an an exponential function (black dashed lines) in order to estimate the characteristic activation time (τ). In the lower pannel, the set of estimated characteristic times was fitted to a Gaussian function, which is our proposed form for the voltage-dependent characteristic time of spHCN gating (Eqs 20 and 21). (TIF)

S3 Fig. Calibration of pH-and voltage-dependent gate of the CatSper channel.
In A and B, experimental data on mouse sperm (circles) are shown, along with model fittings (lines) related to the voltage-and pH i -dependent gate variable in equilibrium,m cs . Model fittings in mouse and sea urchin cases are plotted with dashed and bold lines, respectively. The experimental data of A correspond to the current amplitude of divalent ions (symmetrical Ba 2+ ) produced by taking the holding voltage from 0 mV to −100 mV, under different pH i values, and measured by whole-cell patch-clamp in mouse sperm, (data extracted from Fig. 4c of [72]). Mouse data displayed in A and B were simultaneously fitted to the Eqs 44 and 40. To calibrate sea urchin'sm cs , we took the parameters obtained with mouse data as a starting point and manually adjusted them, constraining that the pH i sensitivity should be within the physiological pH i response observed in sea urchin sperm (area marked in gray). C and D correspond to the G/V (conductance/voltage) curves of mouse and sea urchin CatSper, respectively, using the Eq 39. In C, the parameters reported in [72] were used). In panel D, the resting membrane potential is indicated with a gray dashed line as reference. Taking into account that the S4 segment of CatSper voltage sensor domain has more positive charges in the sea urchin homologue protein that in the mammalian counterparts [12], we envisioned that its voltage sensitivity should be steeper. Thus, a 3-fold decrease was introduced to the voltage sensitivity parameter, s 2 , as an initial guess. This last parameter has not been estimated in sea urchin due to the lack of patch-clamp measurements.  θ). For reference, the original CatSper-only scenario is shown first, i.e. θ = 0, and the subsequent rows correspond to the gradual increase of θ. The effect of this parameter on the channel densities is reported in the coefficient r g , which measures the total conductance ratio CatSper/Ca V , as well as the corresponding g bk value. The greater the value of r g , the greater the predominance of CatSper on the total calcium conductance with respect to Ca V . (TIFF)

S5 Fig. Titration of Ca V or BK channels on the CatSper module.
In A and B, we show numerical solutions for calcium under different values of the weighting parameter θ, with either BK or CaV conductance density set to 0, respectively. As reference, the original scenario of only CatSper, i.e. θ = 0, is shown in the top row, while the subsequent rows correspond to the gradual increase of θ. In A, the reference value of the Ca V conductance density, g cv , is multiplied directly by θ, whereas the reference value of CatSper conductance density is multiplied by (1-θ); The resulting ratio of CatSper/Ca V conductance densities, r g , is shown for each θ value. In B, θ sets the fraction of the reference value of g bk by multiplying the conductance density of BK; unlike scenario A, the conductance density of CatSper is not weighted by θ. For each titration, the corresponding modified g bk value is shown. (TIFF) S1 Algorithm. Two-phase parameter space exploration algorithm.