In the following analysis, all properties have been obtained considering a model with only one pump. Moreover, Table shows that as the cell volume grows, so will be size of the model and the time required to build and check it. We have used a cell volume equal to 10

^{–20} liters for analysis in the following sections. This type of abstraction strategy is common for modeling biological systems as discussed in [

14].

| **Table 2**Model size ranging the cell volume. |

Discovering rare events

Uncommon events can have a significant impact in any system and particularly in biological systems. For example, if a particular combination of reagents can cause a pump to block permanently, it can cause cell death. No matter how unlikely this event is, if it happens the consequences are critical. Traditional analysis methods such as stochastic simulations can miss uncommon or rare events, because they simulate random paths in the evolution of the system, and if the event is rare, it is not likely that it will be simulated in a viable amount of time. PMC, however, can identify these events by looking for them. By stating a property that is true if such an event occurs, PMC can identify the conditions for its occurrence, and as a consequence, uncover hidden but potentially important behaviors in the system.

Our first analysis shows how model checking can be used to identify rare events in the Na,K-pump. Figure presents the potassium concentration in

*M* outside the cell over time for the ODE approach (dashed line and

*y* axis on the right), which uses a deterministic and continuous pump model. The model was built and solved in MATLAB as described in [

15]. The figure also presents the count of potassium ions outside the cell given by a simulation trace (solid line and

*y* axis on the left) of the discrete and stochastic pump model. The simulation trace was obtained using the BIONETGEN tool [

16], which provides an implementation of the direct method of Gillespie.

As can be seen in the ODE approach, the potassium outside the cell is decreasing until around 2 seconds, and its concentration converges on about 0.0018 M. However, this average behavior hides some important traces, as it is shown in the same figure, where the potassium count outside the cell, after the fast decrease until about 2 time units, will oscillate around 12 and might end, i.e. it can eventually reach value 0. However, the probability of this event *potassium outside the cell ends* is extremely low (6.33 × 10^{–3}) during the first 10 seconds. This probability value was determined using the CSL property *P* =? [*F* <= 10 *kOut* = 0]. Notice that there is a significant difference between a rare event and an impossible event. If the event is rare we may decide, based on its occurrence probability, to investigate it further or not. If it is an important event that may cause death it may merit further studies, whereas if it *cannot* happen, we need not worry about it. Properties such as these can help direct future researches and shorten scientific discovery times.

Thus, whereas a deterministic simulation *will* not identify this rare event *potassium outside the cell ends* because it captures only average behaviors, the stochastic simulation approach *may* not capture traces where it happens, depending on the simulation time and on the number of simulated traces.

However PMC can provide stochastic simulation with some hints in this sense. As it lets us know in advance that the rare event happens with probability equal to 6.33 × 10^{–3} in the first 10 seconds, if the stochastic simulation time being considered is 10, in a sample of 1000 traces, for example, about 6 or 7 of them will probably show the rare event.

It is also important to note that in PMC the time is continuous, while in stochastic simulation it is discrete. Hence, if the duration of an event of interest is smaller than the time step being considered in the simulation, it will not be captured, whereas it will be considered in the PMC model. As will be shown, PMC can give some clues for stochastic simulation in order to address these issues. The CSL property shown below (followed by the model checker answer)

ensures that in all model traces the potassium outside the cell, in fact, will eventually end. Additionally, in order to know about the expected time for this event to happen, properties (6) and (7) can be used:

Property (6) means *What is the expected time to the potassium outside the cell being over?*. The model checker answer was 1287 seconds. The reward structure for this property associates reward 1 with each state (see the reward structure ”time” in Sect. *Methods - Sodium-potassium pump specification - PRISM Specification*). On the other hand, property (7) asks *about the probability of potassium outside the cell eventually being over in the first 1287 seconds*, whose answer is 0.63. Thus, we can conclude that during the first 1287 seconds the potassium outside the cell ends in around 63 % of the traces of the system. In this case the expected time is very long and unlikely to be reached before another event occurs. However, the same technique can be used to model reaction time to the presence of a toxin, for example, and the expected time can be crucial to the survival of the cell.

Finally, the following properties are used to reason about the

*maximum and minimum expected time for the system to go from a state where potassium outside the cell is over to a state where this species is not more over*:

There is more than one state where potassium outside the cell is over and, therefore,

**min** and

**max** are used to return the minimum and maximum expected time to reach a state where

*kOut* > 0, ranging over all start states that satisfy

*kOut* = 0. The model checking answers for (8) and (9) properties are, respectively, 111 milliseconds and 14 milliseconds. This minimum expected time for the event duration can be used as a guideline for choosing the time step in stochastic simulation to guarantee that no such events will be lost. Another event that can be easily identified with model checking is if

*after that potassium outside the cell is ended*,

*its amount will eventually return to the initial count KO*. Property (10) is used to verify if this event happens in all traces in the model, whereas property (11) is used to determine the expected time for the event to happen:

where ”kOutOver” is a label to *kOut* = 0. The model checker result for property (10) was true and the expected time computed in property (11) is about 132, 515 seconds.

Thus, two significant events in the system *potassium outside the cell will end and the potassium outside the cell will return to its initial count* will happen in all traces of the pump model and lead us to study about their recurrence in the long term.

Reversibility of the sodium-potassium pump

Due to the fact that there are backwards and forwards transitions for all reactions involved in the Na,K-pump mechanism, as is shown in Fig. , and that many of these transition rates depend on transmembrane substrates, the pump mechanism is automatically reversible, i.e the reactions can be run either forwards or in reverse direction, depending on changes in the amount of substrates. Thus, given the initial concentrations of the substrates, we can consider that the pump performs two main steps. First it runs following the forward reactions and reaches the configuration where substrates reach a maximum or minimum concentration. Then, given these changes in the amount of substrates, the pump returns to the initial configuration through the reverse reactions. Of course, it is possible that the reverse reactions can be followed during the first step, and, similarly, forward reactions can be followed during the second step. In the forward running, Na,K-pump uses one ATP to perform the electrogenic exchange of 2 potassium ions from outside to inside the cell in exchange for 3 sodium ions from inside to outside the cell. In the reverse direction ATP can be produced from ADP and P_{i}.

Without loss of generality, we will study this reversible pump behavior in terms of the potassium amount outside the cell, which will be the species under observation. We can see this pump reversibility as an infinite oscillation between two values, the initial amount of potassium outside the cell,

*KO*, and the final amount of substrates, after the forward running is complete and before the reverse running starts. In this final configuration, the potassium outside the cell ends, i.e. it reaches its minimum value (

*kOut* = 0). This pump reversibility is expressed through CSL property (12), which means

*What is the probability that kOut oscillation between 0 and KO values will never terminate?*where *i* is *KO* and *j* is 0. The result of CSL property (12) is 1, which proves that the events *potassium outside the cell ends* and *potassium outside the cell reaches the initial amount* happen infinitely often, i.e. the automatic pump reversibility property of the pump is true. This pump property can not be seen in the Fig. , because the study of the average behavior of the pump overlooks some aspects of its reversibility. Moreover we are reasoning about the pump infinite behavior, which cannot be achieved through generating and analysis of finite-time trajectories with stochastic simulation. Property (12) can be used to check if the concentration of some species oscillates between any two values *i* and *j*.

Understanding the pump cycle

In this section we present a study of the Na,K-pump mechanism in terms of the rates in the cycle shown in Fig. , in order to understand why the depletion of potassium outside the cell and consequently the pump reversibility can take long periods of time to be completed.

We now introduce some definitions and extensions in the previous PRISM pump model that will be used later. First, we compute the positive or ascending trend [

17] of a transition rate

*r*(

*s*,

*s*_{i}) from a current state

*s*:

where *E*(*s*) is the exit rate of state *s*, i.e., *E*(*s*) = *∑*_{s′S} r(*s*, *s′*), being *S* a finite set of states, and *ξ* is a threshold that indicates a positive trend. We have chosen the value 0.6 for *ξ* in our analysis and, therefore, informally an ascending trend for a transition rate *r*(*s*, *s*_{i}) means that the probability of the system goes from *s* to *s*_{i} (at least 0.6) is greater than goes to any other state *s*_{j} (*j* ≠ *i*), whose transition rate *r*(*s*, *s*_{j}) > 0.

Thus, we add trend formulas to the previous PRISM model for all transition rates using PRISM resources (

*labels* and

*formulas*). The code in

*PRISM Model Extension 2* illustrates the procedure for computing the positive trend label for the transition rate

*r*_{1}.

The rate transition

*r*_{1} is computed by the formula

**rate_r1** and it is different to 0 when the current pump state is

*E*_{1}. ATP and there is enough sodium inside the cell

. In this case, the final value for

*r*_{1} is determined in the same way as described in Sect.

*Methods - Sodium-potassium pump specification - PRISM Specification*. The other rates are computed in similar way than

*r*_{1} and the formula exit rate represents their summation. The probability that

*r*_{1} is taken in the current state is given by formula

**rate_r1_d**, whereas the label

**trend_r1_up** represents if

*r*_{1} really has an ascending trend, i.e. ↑

*r*_{1} is 1. Now, we can use the CSL property (13) to identify the rates that never have a positive trend during the system evolution and, consequently those rates that always have an ascending trend:

In our pump model, [*r*] is used to form the trend label of the backward rates, and *i* ranges from 1 up 6, because there are six reactions in the pump system, considering each direction (forward or backward). Thus, the trend ”trend_r1” represents the trend label for the first pump reaction in the forward direction (see Fig. ), whereas the trend ”trend_rr1” refers to the trend label for the first pump reaction in the backward direction.

The results for property (13) are summarized in Fig. , which shows the trend for all transition rates in the pump cycle presented in Fig. . The arrows leaving the pump states (such as *E*_{1}.ATP) are labelled with the rates for the transition between the current state and the next state, which is given by the direction of the arrow. Associated with each arrow, there is also a sign that indicates if the transition rate has always a positive trend (+) or a negative trend (-), and, finally, if the trend can be negative and positive during the system evolution (+/-).

We can see that the forward rates *r*_{1} and *r*_{2} always have a negative trend, while *r*_{6} always has a positive trend during the system evolution. Moreover, the trends for the forward rates *r*_{3}, *r*_{4} and *r*_{5} can be positive or negative, depending on the changes in the amount of substrates during the pump evolution.

In order to identify the moment when these forward transition rates which don’t have only a positive or negative trend during the system evolution change their trends, we have again extended our PRISM model with the following transition-rewards

In the

*PRISM Model Extension 3*,

**plusKout** is a reward that assigns 1 to each transition from the state K

_{2}.

*E*_{2} to state

*E*_{2} ~ P, which results in the releasing of two potassium ions outside the cell. On the other hand

**minusKout** is a reward that assigns 1 to each transition from the state

*E*_{2} ~ P to the state K

_{2}.

*E*_{2}, which results in the consumption of two potassium ions outside the cell. CSL property (14) determines

*the expected count of potassium outside the cell when the rate r*[

*r*]

_{i} starts to have a positive trend:

Using property (14), we can see that *r*_{3}, *r*_{4} and *r*_{5} (forward rates) start to have a negative trend only when the potassium outside the cell is, respectively, 21, 7 and 7 (the initial amount of potassium outside the cell in our model is 61).

Thus, we can divide the pump operation into three main steps, as is shown in Fig. . Initially (A), despite the general trend to move backwards due to the positive trends of the backward rates *rr*_{1} and *rr*_{6}, once *r*_{3} is taken, the system might complete easily the cycle in the forward direction, because the forward rates *r*_{3}, *r*_{4}, *r*_{5} and *r*_{6} have a positive trend in the most of the time. The backward rate *rr*_{2} needs that the pump goes in the forward direction awhile, increasing the amount of ADP inside the cell, in order to exhibit a positive trend. When potassium outside the cell reaches the value 21, the rates *r*_{3} and *rr*_{2} changes their trends, starting the intermediate step (B). In this step, the pump can still move in forward direction. The last step (C), starts when the potassium outside the cell reaches the value 7, causing changes in the trends of the forward rates *r*_{4} and *r*_{5}. First, rate *r*_{4} no longer has a positive trend, while the negative trend of the backward rate *rr*_{3} is replaced by a positive one. This happens due to the increase of sodium outside the cell, which gives strength to *rr*_{3}, and the decrease of potassium outside the cell, which weakens *r*_{4}. Second, the forward rate *r*_{5} also stops exhibiting a positive trend, whereas the trend of the backward rate *rr*_{4} starts to be ascending. This change is caused by the accumulation of P_{i} inside the cell and the reduction of ATP due to the pump movement in the forward direction. In step (C) there is a low probability, although is not impossible, that the pump continues its operation in the forward direction, given that the only forward rate with positive trend is *r*_{6}, delaying the depletion of potassium outside the cell. In fact, there is a strong general trend for the pump to move backwards, returning to the intermediate step, where the system stays most of the time. Additionally, the pump can move backwards from the intermediate step, returning to the initial configuration. However, this takes long periods of time, given that it is necessary to move against the positive trends of the forward rates *r*_{3}, *r*_{4}, *r*_{5} and *r*_{6} in the initial step.

As shown in the previous sections, the depletion of potassium outside the cell and the pump reversibility are events that can happen in the pump model. However, they can take longs periods of time to be completed. So the study of this section is important to indicate the reasons for this delay. For example, it is possible to see that the first obstacle in the normal operation of the pump is the accumulation of ADP inside the cell which causes the reversion of the *r*_{3} trend. This may indicate a specific aspect of the system that merits further studies. This result may lead to a more precise study because it tells us in detail what has happened (accumulation of ADP inside the cell) and not simply that the pump has reversed its behavior. Results such as these can uncover important hidden behaviors that can speed up further experiments and increase their accuracy.

Validation of the PRISM model

In this section we will show that our PRISM model can produce similar results when compared to the stochastic and deterministic simulations. Property (15) allows us to know the expected amount of potassium outside the cell in time

*T*, which specified in the property:

The label ”kOut” is a reward name defined as shown in *Rewards to PRISM Model, Sect. Methods - Sodium-potassium pump specification - Prism specification*. PRISM supports *experiments*, which is a way of automating multiple instances of model checking. In our case, this is done by ranging the constant *T* from, for example, 0 up to 10, with steps of 0.25. The resulting graph is shown in Fig. (dashed line), which is very similar to deterministic curve shown in Fig. . We also got a similar trajectory using the PRISM tool, Fig. (solid line), which besides verification can also perform stochastic simulation of the model that mimics the Gillespie method. Thus, we can see that PRISM results for the sodium-potassium pump are very close to those obtained using the traditional approaches.