In Silico Determination of the Effect of Multi-Target Drugs on Calcium Dynamics Signaling Network Underlying Sea Urchin Spermatozoa Motility

The motility of spermatozoa of both Lytechinus pictus and Strongylocentrotus purpuratus sea urchin species is modulated by the egg-derived decapeptide speract via an oscillatory [Ca2+]-dependent signaling pathway. Comprehension of this pathway is hence directly related to the understanding of regulated sperm swimming. Niflumic acid (NFA), a nonsteroidal anti-inflammatory drug alters several ion channels. Though unspecific, NFA profoundly affects how sea urchin sperm respond to speract, increasing the [Ca2+]i oscillation period, amplitude, peak and average level values of the responses in immobilized and swimming cells. A previous logical network model we developed for the [Ca2+] dynamics of speract signaling cascade in sea urchin sperm allows integrated dissection of individual and multiple actions of NFA. Among the channels affected by NFA are: hyperpolarization-activated and cyclic nucleotide gated Na+ channels (HCN), [Ca2+]-dependent Cl− channels (CaCC) and [Ca2+]-dependent K+ channels (CaKC), all present in the sea urchin genome. Here, using our model we investigated the effect of blocking in silico HCN and CaCC channels suggested by experiments. Regarding CaKC channels, arguments can be provided for either their blockage or activation by NFA. Our study yielded two scenarios compliant with experimental observations: i) under CaKC inhibition, this [Ca2+]-dependent K+ channel should be different from the Slo1 channel and ii) under activation of the CaKC channel, another [Ca2+] channel not considered previously in the network is required, such as the pH-dependent CatSper channel. Additionally, our findings predict cause-effect relations resulting from a selective inhibition of those channels. Knowledge of these relations may be of consequence for a variety of electrophysiological studies and have an impact on drug related investigations. Our study contributes to a better grasp of the network dynamics and suggests further experimental work.


Introduction
Fertilization is an important process in life. Reproductive success is attained by means of different strategies that increase the probability of gamete encounter. Several species, including sea urchins, produce spermatozoa with swimming patterns regulated by egg secretions. Strongylocentrotus purpuratus and Lytechinus pictus sea urchin spermatozoa swimming is modulated by speract, a decapeptide contained in the outer coating of the egg which diffuses in the sea [1,2]. When these sperm detect speract by means of receptors along the flagellum, an intracellular signaling pathway that regulates fluctuations of the intracellular Ca 2+ concentration ([Ca 2+ ] i ) is triggered. The pathway involves production of cyclic nucleotides and alterations in the ionic permeability of the plasma membrane ( Fig 1A). These biochemical changes are timely translated to enhance the encounter of the sperm with the egg. Sea urchin spermatozoa swimming close to a surface describe circles, a very convenient trait to do imaging [3][4][5]. In the presence of speract gradients a cross-correlation of [Ca 2+ ] oscillations and path curvatures takes place [3,[5][6][7]. Increases in the [Ca 2+ ] i are associated with sharp turning events (high path curvature) that are interspersed with periods of straighter swimming episodes (low path curvature). This swimming pattern is common to a wide variety of organisms with external fertilization [3,4,[8][9][10][11][12][13][14][15][16].
Niflumic acid (NFA), a nonsteroidal anti-inflammatory drug, is able to block or modify several ion channels. Its lack of specificity, usually disadvantageous, turns out to be key to the profound effects it generates on how sea urchin sperm respond to speract. Immobilized S. purpuratus sperm exposed to NFA respond to speract with [Ca 2+ ] i fluctuations that are larger, longer and with increased time intervals between them than control speract responses [17] (Fig 2). Furthermore, in S. purpuratus swimming sperm these alterations on flagellar [Ca 2+ ] i dynamics caused by NFA increase the speractinduced flagellar asymmetry resulting in more pronounced and prolonged sharp turns [4]. In L. pictus swimming sperm NFA inhibits chemotaxis (also in Arbacia Punctulata [15]), without altering their ability to detect the chemoatractant gradient [8].
It is because of the findings mentioned above that it is very important and interesting to unravel how NFA achieves these responses. There are not many tools that efficiently allow analysis of multiple variables regulating a transduction path; certainly a network model of a signaling pathway is one. In this regard using this approach to examine how NFA may achieve its remarkable effects on the speract signaling cascade seems warranted. This Boolean network model has already shown its ability to uncover and predict properties of this signaling cascade. In a previous publication we presented a theory-experiment investigation where the model predicted the possible participation of a CaKC channel and a voltage dependent Ca 2+ channel. Experiments simultaneously recording sperm [Ca 2+ ] i and swimming trajectories, and using specific blockers for these channels gave results consistent with the model predictions [18].
Here we use our model to explore the effect of NFA on the speract signaling cascade of S. purpuratus and L. pictus sperm, always keeping in mind previous experimental results. We focus on Signaling pathway operation diagram, black arrows correspond to activation, red lines to deactivation and yellow arrows can be activating or inhibitory depending on the relative state of the pathway elements being interconnected. Once speract binds to its receptor the several feedback loops are triggered according to the nature of the links involved. The concatenation of these loops leads to oscillatory stages of the whole pathway. The color code identifies corresponding upper and lower part components. Current models propose that the binding of speract to its receptor promotes the synthesis of cGMP that activates K z selective and cyclic nucleotide-gated channels (KCNG) leading to membrane potential (V) hyperpolarization [3,4,[7][8][9][10][11]18]. This V change first induces an intracellular pH increase via a Na z /H z exchanger (NHE) activation, [18,51,52], stimulates hyperpolarization-activated and cyclic nucleotide-gated channels (HCN) [53][54][55], removes the inactivation of voltagegated Ca 2z channels HVA and LVA [18,56] (CaV), and facilitates Ca 2z extrusion by Na z /Ca 2z exchangers (NCE) [51,52]. The opening of HCN and the influx of Na z contribute to V depolarization, and concomitant increases in ½Ca 2z i and ½Na z i further depolarize V. It has been proposed that the ½Ca 2z i increases could lead to the opening of Ca 2z -regulated Cl { channels (CaCC) and/or Ca 2z -regulated K z channels (CaKC), which would then contribute to hyperpolarize the V again, removing inactivation from CaV channels and opening HCN channels, [3,4,18]. It is thought that this series of events is then cyclically repeated generating a sequence of V-dependent turns. B) Network model of the signaling pathway. The network can be envisaged as a circuit where each node represents an element of the pathway and links, either in the form of arrows or lines, correspond to connections determined in the bottom part of (A). The activating or inhibitory nature of the yellow lines depends on the value of voltage (V). Yellow nodes represent binary nodes (0,1), and the four brown nodes are ternary nodes that can take values 0, 1 and 2. Changes in the node states are determined by the connected nodes by means of a regulatory function (or truth table). As an illustration we present the case of the cGMP shown at the bottom left of (B). The first 3 columns in this table contain all the possible activation states of the cGMP regulators: GC, which is an activator; PDE, an inhibitor and cGMP (cGMP is a self-regulator); the fourth column shows the values for the function that correspond to each combination of the regulators. Additional nomenclature note: Speract receptor (SR); guanylate cyclase (GC); unknown Ca 2z channels sensitive to cAMP (cAMPCC); Ca 2z pump ( four individual-cell quantities measured previously by Wood et al., in [4], namely, average [Ca 2z ] i level, amplitude, peak and frequency. We use them to compare the model predictions between oscillations in the wild type (WT) untreated network and in the NFA treated case, considering the two scenarios, depending on wheter NFA inhibits or activates CaKC channels.

Methods and Models
The logical network The logical signaling network corresponding to the SASP, first introduced in [18], consists of nodes interlinked according to Fig 1A, representing the principal components involved in the signaling cascade: ion channel activities, intracellular ion and molecular concentrations and the membrane potential amongst others. The network is shown in Fig 1B and the nomenclature is explained in its figure caption. To analyze the dynamics of the network, we implemented a discrete formulation that is a generalization of the Boolean approach and that has proven to be revealing for the gene regulation dynamics of many systems [38][39][40][41][42], as well as other cell signaling networks [43]. In this approach, the dynamical state of the network consists of a set of N discrete variables fs 1 , s 2 , . . . , s n g, each representing the state of a node. For this particular network, most of the variables take on two values, 0 and 1, depending on whether the corresponding element is absent or present, closed or open, inactive or active, etc. However, an accurate description of the dynamical processes in the network required four nodes to be represented by three-state variables: the membrane potential (hyperpolarized 0, resting 1, and depolarized 2); the low and high threshold voltage-gated Ca 2z channels (inactive 0, closed 1, and open 2); and the intracellular [Ca 2+ ] concentration ½Ca 2z i (basal 0, tonic 1 and supratonic 2). The state of each node s n is determined by its set of regulators (which are some other nodes that also belong to the network). Let us denote as s n1 , s n2 , . . . , s n k the k regulators of s n . Then, at each time step the value of s n is given by where F n is a regulatory function constructed by taking into account the activating/inhibiting nature of the regulators. Each node has its own regulatory function. For the construction of these regulatory functions, we have made use of extensive biological knowledge, mainly of an electrophysiological nature, available to us in the literature and in our own laboratory. An illustration of such regulatory functions is given in Fig 1 for the cGMP node. For the case of [Ca 2+ ], which is the main concern in our study, we have that it is regulated by 6 nodes (see origin of incoming arrows and red lines of dCa in Fig 1B). LVA, HVA and cAMP-dependent Ca 2+ channels are activators, i. e., their activation at time t, will favor and increase in [Ca 2+ ] at the subsequent time t + 1. There are also two inhibitory nodes: the Ca 2+ pump and the Na + /Ca 2+ exchanger which when activated at time t will favor a decrease in [Ca 2+ ] at time t + 1. Finally the sixth node is the [Ca 2+ ] node itself that acts as an inactivator, hence favoring a decrease in its value at the following time step. The balance among all the above inputs determined by physiological considerations is expressed in the form of a truth table with 432 rows, which constitutes its regulatory function. In the supplementary material S1 the regulatory functions for all the networks nodes are shown explicitly. With this model we can therefore observe in silico the effect of altering certain elements relevant to the pathway. In this paper we consider the case of NFA-sensitive channels: HCN, CaKC or CaCC. In order to test the effect of NFA on the network evolution, we proceeded with the deletion or activation of these channels according to the two scenarios mentioned above, individually, by pairs (HCN

Steady state analysis
For this type of network (finite number of nodes which take a finite number of values), all initial conditions lead to a periodic behavior where the network configuration is replicated after a certain number of steps. The time required to reach this condition is known as the transient time; it is important to mention that, independently of the initial condition, it is shorter than 45 steps. The attractor is reached after 22 iterations and the number of iterations between the repeated configurations is the period. These periodic solutions are the attractors of the network dynamics. For the wild type (WT) network (with all nodes present), 88.9% of all speract activated initial conditions lead to a period-4 attractor, while the remaining 11.1% converge to a period-8 attractor [18]. The total amount of initial conditions which reach a particular attractor is called the basin of attraction. In order to understand the differences between the averaged [Ca 2+ ] concentration overline½Ca 2z i steady-state behavior of NFA-treated network and the WT, we calculate modifications of the network attractors after alteration of all combinations of NFA-sensitive channels, paying attention to the ensuing number of attractors and their associated periodicity. We make the comparisons between the characterization of the periodic behaviors of ½Ca 2z i and the period of the network attractors. This is done by means of a Fourier spectrum analysis of ½Ca 2z i from which we can determine differences in the temporal behavior of WT and NFA-treated network dynamics. We should point out that we have explored the effect of small perturbations on the regulatory functions, such as modifications in the outcome of any row in the truth tables of all the nodes. This holds particularly well for [Ca 2z ] and is a measure of the high robustness of the dynamics we are implementing and hence of the reliability of our results.

Individual Channel Results
In the first column of Fig 3,  ] node network dynamics period-4 attractor together wtih the second harmonic of the period-8 attractor. We shall say that the ½Ca 2z i time series has a richer temporal behavior if its Fourier decomposition is more complex. Blockage of the HCN channel in our network model produces a more elaborate temporal behavior, which is registered in the Fourier power spectrum by the appearance of multiple lines that correspond to different modes in the Fourier decomposition. An opposite effect is observed with the blockage of the CaCC channel, where the Fourier decomposition consists of only one red line at the frequency value of 1/4 that corresponds to a unique period-4 attractor of the [Ca 2z ] dynamics. For the case of CaKC blockage the Fourier spectrum remains unchanged with respect to the WT; this is consistent with the persistence in the blocked network of the WT attractors. Finally, activation of the CaKC channel produces 3 peaks at 1/7, 2/7 and 3/7 frequency values which are related to a period-7 attractor obtained by the network dynamics.

Double Channel Results
With regard to the temporal behavior, a similar analysis to the one performed for Fig 3, reaffirms the trends found from the individual channel alterations (Fig 4), namely: i) HCN blockage leads to a more complex Fourier Decomposition, consistent with a more elaborate behavior of the ½Ca 2z i series shown on the left. ii) CaCC blockage reduces the Fourier components leading to a simpler periodic behavior. iii) CaKC activation enhances the effect of the concomitant deletion of either HCN or CaCC.

Triple Channel Results
With the joint alteration of the HCN-CaCC-CaKC channels the differences of the effect of blocking or activating the CaKC channel can be appreciated in Figs. 5 and 6 respectively. Fig 5A is a segment of the ½Ca 2z i time series once steady-state conditions have been reached. Fig 5B is the Fourier spectrum calculated over 1000 steady-state time steps. Besides the coincidence of the WT period 4 and 8 peaks with their harmonics as before, a new feature that arises in the dynamics of the blocked system is the occurrence of a period-9 peak and one of its harmonics. This is an indication of a richer temporal behavior. If we perform a running average of the ½Ca 2z i series over a window of 4 steps, the smaller fluctuations are dampened and the behavior of an envelope at larger scales becomes evident. This is shown in Fig 5C where a recurrent module of 72 time steps comes to light. Note that 72 is the minimum common multiple of the period 8 and 9 Fourier components, this is a manifestation of a period locking phenomena at the level of the ½Ca 2z i network attractors. If the running average is performed now over an 8 step window, then the WT ½Ca 2z i fluctuations are completely leveled out, as well as the period-8 component of the blocked case network, allowing for the period-9 component to surface (see Fig 5D). When activation of CaKC channel is considered, we encounter the behaviors shown in  Fig 6B. When a 4 step average is performed, the WT is strongly smoothened and a period 8 structure emerges while the structure of NFA treated case remains similar (see Fig 6C). In Fig 6D we show that though the WT is completely flatten to its over-all average, under an 8 step running average, when the effect of NFA is taken into account a module 3 structure is preserved. For the converse case of averaging over 3 steps, it is the treated case that gets flatten while the 8 module consisting of two alternating 4 modules structure of the WT survives.

Overall Picture
Given our in silico results, we can conclude from our data, the following: Blockage of the HCN channel affects mainly the temporal behavior of the [Ca 2z ] dynamics. This became evident in figures 3 and 4, where changes in the temporal structure of the Ca 2z curve as well as the number of Fourier components can be observed. Another visible effect is the reduction of the ½Ca 2z i amplitude in all cases where HCN participates. When CaCC is blocked, the main effects are the increase of the amplitude in the ½Ca 2z i curve and a recovery of the regularity. Fourier components are reduced in all cases in which CaCC is altered. Regarding scenario 1, in which CaKC channel is inhibited, the principal effect is an elevation in the ½Ca 2z i mean and maximum peak. In scenario 2, after activating the CaKC channel, the overall effect is quite similar to the one observed by HCN elimination: reduction in ½Ca 2z i amplitude, peak and mean but, an increase in Fourier modes. Though all of these effects are reproduced under combinations of the above channel alterations, it is possible to establish the dominance of one channel alteration over another. For example, with regard to the fluctuation amplitude we have that the effect of deleting HCN overrides CaCC blockage, however, there is a recovery of regularity, situation in which the effect on temporal behavior of CaCC blockage overrides HCN elimination (Fig 4A). The same temporal observation is obtained by looking at the Fourier power spectrum for the case of the joint blockage of HCN and activation of CaKC. Both alterations result in an increase of the Fourier modes (Fig 4K) not only with regard to the WT but also to their action taken separately (Figs. 3C and 3L). This last result has a physiological explanation, because the effect of blockage of a channel which allows a cation entrance (HCN) or the activation of a channel which allows the cation extrusion (CaKC) will have a similar result in the regulation of membrane potential: shorter depolarizations.
Finally, the effect of altering the three channels together is analyzed. Here, we encounter the two scenarios mentioned previously, depending on whether CaKC channel is inhibited (Scenario 1, Fig 5) or activated (Scenario 2, Fig 6).
Overall our main findings are the following: With respect to scenario 1, our analysis shows that only when NFA blocks simultaneously the three channels under consideration, we recover the experimental observations of [4] with regard to WT conditions for the model ½Ca 2z i time steps, namely: a. The mean taken over the time increases (Fig 5A, C and D). b. The amplitude grows (Fig 5D). c. The peak values are higher (Figs 5A, C and D). d. The time interval between oscillations increases (Fig 5B).
In scenario 2 (Fig 6), the effect observed is a decrease in amplitude, period, and maximum peak, a trend opposite to the experiments. This behavior suggests the need of another [Ca 2z ] channel in the signaling pathway. An integrated presentation of our findings is condensed in Table 1. The red arrows correspond to those cases in which CaKC channel is activated (scenario 2). Table 1 clearly shows that within our modeling, experimental results from [4] can only be recovered under the joint NFA deletion of the three HCN, CaCC and CaKC channels. The first three columns of Table 1, state the cause-effect connections of individual channel deletions: the blockage of the CaCC channel increases ½Ca 2z i amplitude fluctuation, the blockage of the CaKC increases the ½Ca 2z i peak and mean values, and the blockage of the HCN channel confers a more elaborate temporal behavior.

Summary
We developed a novel analysis of the effect of the multi-target drug niflumic acid, based on a logical network model of the speract-activated signaling pathway. We selectively modified the functionality of three ionic channels which are sensitive to NFA:   [8,44]. Although NFA is known to affect these three channels, we addressed the question of whether or not a concerted alteration takes place in the sea urchin, taking into account that most pharmacological results have been derived in mammalian channels [19,21,23,24]. If CaKC channel is inhibited, we concluded that the experimental observations of [4] are retrieved only under a concomitant action of NFA on all three channels. The degree of synergy of this action remains to be elucidated. If CaKC channel is activated by NFA, the need of an additional [Ca 2z ] channel comes to light. Furthermore, along our investigation, indications of the relative importance of individual channel deletions surfaced. Our model predictions call for the implementation of experimental protocols for their corroboration.
Knowledge of this type of cause-effect relations may be helpful for the understanding of the actions of diverse molecular compounds and might contribute to elucidate operating mechanisms of drugs.

Physiological Context
The physiological effect that NFA may exert on the Ca 2zactivated K z channel CaKC in the speract-activated signaling pathway (SASP), either activating or inhibiting, can be explained in both cases (Fig 7): i) When CaKC is blocked, changes in membrane potential are attained, the depolarization is more pronounced owing to the reduction of ½K z efflux (Fig 7A, violet line), hence, the voltage-dependent ½Ca 2z channels (CaVs) will remain open for a longer period; consequently, CaKC blockage (scenario 1) produces a bigger [Ca 2z ] curve ( Fig 7B, violet line). Since there is pharmacological evidence of the presence of Slo1-type Ca 2z -activated K z channel in the SASP [18], and it is known that the Slo1 channel is activated by NFA [45,46], the above line of reasoning could hold if we consider the participation of a different CaCKC such as SK Ca or another K z channel in the pathway. Given the results of the last column of table 1, these findings support this alternative, hinted by the experimental discrepancy mentioned in the introduction. ii) In the case of scenario 2, after the first KCNG-dependent hyperpolarization, the activation of the CaKC will produce subsequent higher hyperpolarization values since this CaKC channel remains open. This in turn causes shorter lived depolarizations of a smaller magnitude, situation that leads to a premature closing of CaVs with the consequent decrease in ½Ca 2z i . Under the assumption that NFA activates the CaKC channel, given that NFA is a non-specific drug, the inconsistency of the results between the ½Ca 2z i network dynamics with a CaKC channels activated in silico and experiments, favors the consideration of other [Ca 2z ] channels in the SASP sensitive to NFA, so far not present in our model. Suggestions along these lines are met by the sperm-specific CatSper channel, present in the sea urchin genome [47][48][49], which has recently been found to be a polymodal chemosensor in human sperm [50]. The pertinence and validity of this last suggestion, is a subject for further study. . If the effect of the drug is an inhibition of CaKC channels, this would cause a bigger depolarization, because a decrease in the ½K z i efflux. This is shown in the purple curve, notice that the amplitude is bigger than in the WT case. If the effect on CaKC is an activation instead, a bigger hyperpolarization is produced, due to the increase in the efflux of ½K z i . Depolarization takes less time and is smalles than in the WT case for the same reason (blue curve). For (B), the WT [Ca 2z ] dynamics is again depicted in black. Inhibition CaKC channels (purple curve), reduces the extrusion of ½K z i diminishing hyperpolarization hence delayning the closure of CaV channels. This will cause the elevation of ½Ca 2z i as well as the time between [Ca 2z ] peaks (the period). Overall, the activation of CaKC channels (blue curve), produces a shorter [Ca 2z ] oscillations due to the faster hyperpolarization, which in turn closes the CaV channels avoiding a big elevation of ½Ca 2z i . doi:10.1371/journal.pone.0104451.g007