Mathematical model reveals that heterogeneity in the number of ion transporters regulates the fraction of mouse sperm capacitation

Capacitation is a complex maturation process mammalian sperm must undergo in the female genital tract to be able to fertilize an egg. This process involves, amongst others, physiological changes in flagellar beating pattern, membrane potential, intracellular ion concentrations and protein phosphorylation. Typically, in a capacitation medium, only a fraction of sperm achieve this state. The cause for this heterogeneous response is still not well understood and remains an open question. Here, one of our principal results is to develop a discrete regulatory network, with mostly deterministic dynamics in conjunction with some stochastic elements, for the main biochemical and biophysical processes involved in the early events of capacitation. The model criterion for capacitation requires the convergence of specific levels of a select set of nodes. Besides reproducing several experimental results and providing some insight on the network interrelations, the main contribution of the model is the suggestion that the degree of variability in the total amount and individual number of ion transporters among spermatozoa regulates the fraction of capacitated spermatozoa. This conclusion is consistent with recently reported experimental results. Based on this mathematical analysis, experimental clues are proposed for the control of capacitation levels. Furthermore, cooperative and interference traits that become apparent in the modelling among some components also call for future theoretical and experimental studies.

[4], Figure 6, with some modifications and additions. The gray bands represent the flagellum membrane. The gray semicircle represents the internal calcium reservoir. Inserted inside the light blue background (intracellular medium), we can see from left to right circles representing ion concentrations (green) and membrane potential (V, red). From left to right, inserted inside the upper gray bands: cholesterol acceptor (ChAcc), sodium/calcium exchanger (NCX), calcium pump (PMCA), voltage/pH dependent calcium channel (CatSper), voltage/pH dependent potassium channel (Slo3), electroneutral but voltage dependent sodium/proton exchanger (sNHE). From up to down, inserted inside the gray semicircle (Redundant Nuclear Envelope, RNE): calcium dependent IP 3 Receptor calcium channel (IP 3 R), calcium pump from the redundant nuclear envelope (SPCA). From left to right, inserted inside the inferior gray bands: electrogenic sodium/bicarbonate co-transporter (NBC), electroneutral chloride/bicarbonate exchanger (SLCAE), electrogenic chloride/bicarbonate exchanger (SLC26), chloride channel (CFTR/TMEM16A) and sodium channel (ENaC/TRP), soluble adenylate cyclase (sAC), (PDE) phosphodiesterase, (cAMP) cyclic adenosine monophosphate, (PKA) protein kinase A, (IP3) inositol triphosphate, Cap is a marker that represents the beginning of late capacitation. Due to scarcity of experimental information on chloride channels and sodium channels relevant to murine sperm capacitation, we modeled CFTR and ENaC only, which are the channels with more available experimental support. Blue sharp arrows correspond to positive regulation; red flat arrows, to negative regulation; green triangle arrows, to an increase or a decrease of ion concentration or membrane potential during capacitation; and black triangle arrows, to ion fluxes direction.
The Boolean model can be generalized to multi-state logical models, in which each node 162 σi is allowed to take m states instead of only two: whereẽi is the set of values σi can take, ei j ∈ Z. 164 Most of the nodes in the capacitation network have two dynamical states. In the case of 165 ion channels, the dynamical states portray: open (state 1) and closed (state 0). In the case . Gray rectangles: auxiliary nodes for concentration and membrane potential recovery (see Sec. 4.2.6). Light pink triangle: membrane potential. Black arrows: positive regulation. Red arrows: negative regulation. Dashed arrows: regulation related to fluxes. Gray arrows: negligible regulation due to our coarse-grained discrete modelling approximation; in the case of the joint node regulators, gray arrows means that these regulators are not included due to the relative scarcity of experimental single-cell measurements. NaC incorporates the sodium channels like ENaC/TRP. ClC incorporates the chloride channels like CFTR/TMEM16A. For more information about auxiliary node KSper flux, see Sec. 4.2.3 they have three dynamical states that represent the current direction: inward (state 1), null 170 (state 0) and outward (state −1). In the case of membrane potential, the dynamical states 171 portray: depolarized (state 1), resting (state 0) and hyperpolarized (state −1).

172
In general, there are many ways reported in the literature for building regulatory functions, wi j σi j (t) , wherein the sgn function discretizes the sum of regulators σi j weighted by wi j , according to: where x is a dummy variable. As in Sec. 2.1.2, we can generalize the neural network scheme 209 to a multitemporal variant scheme as follows: where the nodes σi now have a neuron-like regulatory function in which each regulator σi j    j will or will not contribute in the updating ofĨ r i at time t + 1, depending on the random 246 variable b t ν j , which can be either 0 or 1. Note that every spermatozoon r of the population 247 has its own set of regulatory functions {F r i } for the net ion fluxesĨ r i which can differ from 248 spermatozoon to spermatozoon. Also that regulator nodes ∆µi j (t) represent the product of 249 the emf associated with ion transporter j and the respective gating variable, whereas weights 250 wi j would be the analog of effective conductances. The product of these terms gives an ion 251 current expression based on Ohm's law. Regulator nodes ∆µi j (t) have a state repertoire 252 defined by ∆µi j (t) ∈ {−1, 0, 1}, where −1 corresponds to an emf that generates outward ion 253 flux, 0 corresponds to null force, and 1 corresponds to an emf favoring ion influx. In Figure 3 254 we can see how a change in ion concentration xi occurs through the net ion fluxĨ r Kirchhoff's law. Hence: where zj is the charge sign of the reference ion in the current carried by transporter j. The 267 latter quantity ensures that ion fluxes that are capable of modifying the membrane potential 268 follow the typical convention used in electrophysiology, i. e. negative ion currents correspond where si allocates the proper sign for the selected variable i. Summing up, the joint node tells 302 us if, on average, the four selected variables were in their respective capacitation-associated 303 states for a large enough fraction of time θ in an observation window W before the time t, i. e.

304
at local temporal scale. 305 We will say that a single cell is early capacitated if the steady state time-average of σ r J is 306 above another threshold θc. Hence, we propose a capacitation criterion in terms of the binary 307 classifier Θc (Heaviside function), defined by: where σ r cap is the effective output of our capacitation model, which is a discretization of the In the first column, only cholesterol acceptor (ChAcc) was added. In the second column, only bicarbonate was added. In the last column, both stimuli were added to the extracellular medium. For each box, the upper right corner corresponds to changes observed in time series averaged over a 2 × 10 3 simulated sperm population (e. g. Figure  S2). Symbols −, 0 and + refer to a decrease, no apparent change and an increase, respectively. In addition, lower corners contain references for that specific measurement, or NE if there is no experiment reported in the literature. Green indicates agreement between simulations and available experiments, whereas white indicates no applicable comparison due to lack of empirical evidence. The superscript * refers to transient behavior.

LOF-network variant CatSper
Slo3 Cl --channel * Na + -channel ** PKA   Table 1, HCO e +ChAcc), between sperm with a given blocked node and the reference wild-type (∆ lof wt = σ lof i − σ wt i ). The lower-left corner refers to experimental evidence changes, accompanied by references (if available), and the upper-right corner refers to changes corresponding to time series averaged over a 2 × 10 3 simulated sperm population with a variability of D = 0.25 in the distribution of ion transporters. Entries have the following meanings: NE accounts for cases where no experiment is reported in the literature for the particular observable under the specified node blockage; −, 0 and + correspond to a decrease, no apparent change and an increase with respect to the capacitation response in wild-type control. Colors refer to the level of agreement between experiment and simulations: full coincidence (green), insufficient sensitivity (yellow), and no applicable comparison (white). We consider insufficient sensitivity a form of partial agreement wherein the empirical evidence does show a change after a given loss-of-function perturbation, but the model does not display enough sensitivity (i. e. it remains in resting state) to reflect the expected change. *The most likely candidates for Cl --channels would be CFTR and/or TMEM16A. **The most likely candidates for Na + -channels would be ENaC and/or TRP.  The point we want to highlight is the recurrence and robustness of this phenomenology. 420 We have uncovered a procedure dependent on the degree of variability in the number of ion 421 transporters, that leads to capacitation and that can regulate its fraction in the population.

422
Experimental corroboration of this model scenario is essential. Furthermore, the understand- CatSper, PKA) and those directly modifying electrochemical variables related to capacitation 455 response (Na and Cl channels). In addition to the previously presented LOF scheme (Table   456 2), we explore another perturbation scheme: gain-of-function mutants (GOF), in which the 457 node of interest is fixed to the most active state. After a capacitating stimulus, the resulting 458 capacitation fraction in each mutant is depicted with a colored box. As reference, a population

465
In the triangle of GOF-only network variants, the predominant effect is an increase of the 466 capacitation levels, except for CatSper GOF , ClC GOF and NaC GOF combinations.

467
The mixed LOF and GOF scenarios display a more complex pattern than the two above-468 mentioned types of mutants. More specifically, this kind of scenarios allow us to assess pos-469 sible recovery or interference between nodes. For example, the model predicts that the over- the WT capacitation levels. 474 We remark that since the over stimulation of either Slo3 or PKA can maximize the capac-  On the model building 491 We construct a discrete, mostly deterministic regulatory network for early capacitation. The The definition of a capacitated "state" within our modeling scheme, is crucial for our study.

505
Keeping this in mind, we proposed that the convergence of certain levels of short-range tempo-506 ral averages of the set of selected nodes: PKA, Cai , Nai , pHi and V provides a good working 507 definition. We are aware that as capacitation is a broad complex process, recall that here we 508 have only addressed the early stages.    While emphasis is placed on the importance of capacitation regulation due to the variability 622 in the number of ion transporters, it is also worth to highlight that the model is able to 623 incorporate changes to multiple nodes within the system and to assess whether these changes 624 result in increases in capacitation fraction or decreases, thus allowing for a better exploration 625 of oversized additive or multiplicative effects of components versus cancellation.  In this subsection, we describe the building methods for regulatory functions of: ion concen-630 trations, ion transport, emf nodes, auxiliary nodes and early phosphorylation, among others. 631 We give a table with all initial conditions used in our simulations. Also, particular cases of ion flux nodes with modifiers are explained.

633
Due to the discrete formalism, we employ a discretization of Ohm's law, Nernst potential   Table 3 shows the discrete numerical values used in all our simulation as initial conditions.
where PV is the opening probability of CatSper channels by influence of membrane potential     are both in state 1, Slo3 can take the state 1 even in resting potential. If any of the above 687 mentioned conditions is not fulfilled, then the conductance will return to its basal value 0.
where PNaC is the Na-channel gating node output value, NNaC is the number of Na channels  apparent contradiction, we decided to model this current as a purely K + -current from possible 725 additional K + channels downstream of Slo3 activity and with higher selectivity. In order to 726 obtain some notion of magnitude and direction for this current, we considered a potassium 727 flux according to the K + Nernst potential and Ohm equation as follows: 728 where P Slo3 is the Slo3 gating node output value, N Slo3 is the number of Slo3 channels and [79], we took Ki as a constant during all the simulations. Note that Ki is not shown in Figure   737 2 and Table 4 because those only include nodes that vary with time.
where Hi and He (intra-and extracellular proton concentrations) are calculated as 10 −pH i 748 and 10 −pH e , respectively. We assigned state −1 to the sodium outward flux coupled with to the regulator gradients and membrane potential: We assigned the state −1 to outward flux, state 0 for null flux and state 1 for influx 759 of sodium and bicarbonate. The sodium/bicarbonate stoichiometry per individual transport 760 event by this cotransporter, nNBC, is reported in Table 6.

761
Electroneutral chloride/bicarbonate exchanger. In order to obtain the ion flux direction of Cl − /HCO − 3 electroneutral exchanger (SLCAE), we used the following energetic 763 consideration for determining the most probable state of SLCAE according to the regulator 764 gradients: We assigned the state −1 to influx of chloride and efflux of bicarbonate, state 0 for null We assigned state −1 to influx of chloride coupled with extrusion of bicarbonate, state 0 772 for null flux, state 1 for extrusion of chloride coupled with influx of bicarbonate. The chloride 773 to bicarbonate stoichiometry of this exchanger, nSLC26, is reported in Table 6. there is adenosine monophosphate in cytosol (state 1).

780
Cyclic Adenosine monophosphate (cAMP) has three regulators: itself, sAC and PDE. The 781 cAMP will remain in its previous state if both sAC and PDE remain inactive. There will be 782 production of cAMP (state 1) as long as sAC is activated (state 1). There will be a decrease 783 of cAMP (state 0) only when PDE is activated (state 1). There will be a decrease of cAMP 784 (state 0) if PDE and sAC are both activated (state 1).

785
Protein kinase A (PKA) will remain activate (state 1) while there is cAMP in the cytosol 786 (state 1) otherwise PKA will remain inactivated (state 0). When the recovery nodes of (∆µpHR, ∆µNaR and ∆µHCO3R) turn on, they set their respective 789 ion concentration (pHi , Nai and HCO3i) from the high (state 1) to their resting levels.

790
The leak current node (IL) is a nonspecific counterbalance mechanism to voltage changes,  The inositol trisphosphate (IP3) node will increase its concentration in cytosol (state 1), pro-795 vided that both Cai and cAMP remain high (state 1). We discretized the continuous model 796 proposed by [80], in which it is argued that IP3 is synthesized near the midpiece through a 797 PLC isoform and can be regulated by calcium and cAMP.

798
IP3 sensitive calcium channels from calcium reservoirs (IP3R) node will promote a calcium 799 current into the cytosol (state 1) as long as IP3 remains high (state 1) and Cai remains low

Thresholds of membrane potential node
Eq. 11 was deduced from a discrete version of Hodgkin-Huxley equation: where Cm is the flagellum capacitance of the sperm andĨj(t) is the flux of ion transporter j 810 at time t. In general, the weighted sum in Eq. 11 can take values from the interval x ∈ (lo, lm) 811 which can be partitioned into pj = (lj−1, lj) and discretized as follows: where x is a dummy variable, lj ∈ (lo, lm), j ∈ (1, ..., m) with m number of states that 813 membrane potential can take, lo and lm are the minimum and maximum values the weighted 814 sum can take respectively. In particular, for the case of the membrane potential in Eq. 11, 815 the function Θ discretizes ad hoc the output of the weighted sum as follows:   Table 6), and we define wi j as follows: where, for a given ion channel j, Pj is its gating probability, wi j ∈ R and quantities with    In order to introduce variability in the amount of ion transporters present in the sperm popula-   Chloride fluxes. We determine the chloride flux integrator nodeĨ Cl regulation with: where the chloride transporter modifier functions P t and Q t are defined as follows: wI Cl SLC26 G r SLC26 (1, D)∆µSLC26(t) where the bicarbonate transporter modifiers P t and Q t are defined in the same way of chloride.

900
Sodium fluxes. We determine the sodium flux integrator nodeĨNa regulation: wI Na NaR G r NaR (1, D)∆µNaR(t) where the bicarbonate transporter modifier Q t is defined in the same way of chloride.    Figure Figure S3: Averaged time series with variability, under an external stimulation of the addition of cholesterol acceptors and higher bicarbonate levels in the extracellular medium, of a selected set of variables in a WT network as reference (left column), and in a CatSper LOF network variant (right column). Under CatSper LOF , membrane potential V hyperpolarizes, Ca i goes to basal levels, Na i decreases, whereas Cl i , pH i and PKA activity increase. Population size is 2 × 10 3 sperm, the standard deviation used in introducing variability on ion transporter weights is D = 0.25.   Figure S4: Comparison of the averaged time series with variability, under an external stimulation of the addition of cholesterol acceptors and higher bicarbonate levels in the extracellular medium, of select nodes under WT conditions (first column) and loss of function (LOF) of ClC (second column). Notice that under LOF membrane potential V hyperpolarizes, Na i goes to basal levels, Ca i , Cl i , pH i and PKA activity increase. Population size is 2 × 10 3 sperm, the standard deviation used for variability on ion transporter weights is D = 0.25.   Figure S5: Averaged time series with variability, under an external stimulation of the addition of cholesterol acceptors and higher bicarbonate levels in the extracellular medium, of select set of variables under WT conditions (first column) and NaC LOF (second column). Membrane potential V hyperpolarizes, Na i decreases, Ca i , Cl i , pH i and PKA activity increase. Population size is 2 × 10 3 sperm, the standard deviation used in introducing variability on ion transporter weights is D = 0.25.   Figure S6: Averaged time series of select set of variables with variability, under an external stimulation of the addition of cholesterol acceptors and higher bicarbonate levels in the extracellular medium, under WT (first column) conditions and Slo3 LOF (second column). Membrane potential V depolarizes, Na i decreases, Ca i , Cl i , pH i and PKA activity increases. Population size is 2 × 10 3 sperm, the standard deviation used in introducing variability on ion transporter weights is D = 0.25.   Figure S7: Averaged time series with variability, under an external stimulation of the addition of cholesterol acceptors and higher bicarbonate levels in the extracellular medium, of select set of variables under WT conditions (first column) and PKA LOF (second column). Membrane potential V depolarizes, Na i , Ca i , Cl i and pH i increase, PKA activity goes to zero. Population size is 2 × 10 3 sperm, the standard deviation used for introducing variability on ion transporter weights is D = 0.25.  Figure S8: Capacitated fraction (CF) in sperm populations for fixed joint threshold θ = 0.175, varying three parameters; observation window size (W , in iteration units), standard deviation of Gaussian variability (D) and activity threshold of joint node (θ c ). In all surfaces, each data point comes from 2 × 10 3 -sperm sized populations with 2 × 10 4 iterations long simulations (5 × 10 3 iterations at resting state and 1.5 × 10 4 post stimulus iterations). The classification of capacitated vs. non-capacitated sperm is applied on the last 10 4 time iterations of each individual sperm. The color bar represents the level of capacitation fraction. Note that the white zone corresponds to a capacitation level of 30 ± 5%, which is close to the levels typically observed for in vitro capacitation in wild-type sperm.  Figure S9: Capacitation fraction (CF) surfaces in sperm populations, for fixed joint threshold θ = 0.250, determined as in above Figure S8 .