# Limit Cycle Oscillations in Pacemaker Cells.

## I Introduction

In elementary electrostatics it is well known that the relation between the voltage and the charge of a capacitor is

$$q=Cv,$$
(1)   
where $`v`$ is voltage, $`C`$ is capasitance and $`q`$ is charge. Differentiating this equation with respect to time we obtain

$$\frac{dq}{dt}=C\frac{dv}{dt},$$
(2)   
where $`\frac{dq}{dt}i`$ is a current, and the sign of $`v`$ is a matter of convention. In physiology this second relation (2) is being used to describe how the membrane potential ($`v`$) is changed when ions move across the cell membrane. Unfortunately we observe that this equation also is being used in some models where one in addition keeps track of the charge ($`q`$) concentrations inside and outside the cell . In those models voltage and charge are believed to be independent dynamic variables: first one determines the voltage by integrating the membrane currents, then one determines the charge by integrating the same membrane currents.

The purpose of this article is to point out that integrating the membrane currents once is enough. Voltage and charge cannot simultaneously be independent dynamical variables in a model, simply because of (1).

In order to visualize the drawbacks of treating voltage and charge as independent variables, we explore numerically the nonlinear dynamics of two different models describing the pacemaker activity of the rabbit sinoatrial node. The procedure is as follows:

1. we integrate numerically the equations of motion for a sufficiently long time to detect a steady state,
2. we change the initial conditions and repeat 1.

First we display results from the Wilders et al. model , a model that treats voltage and charge as independent variables. In that model it is thus possible to select an initial voltage and an initial charge independently. The dynamics of that model seems peculiar. An infinity of limit cycles is observed: each time we select new initial conditions a new limit cycle, corresponding to a new value of the constant of motion $`qCv`$, is found. This hampers the usefulness of the model.

Second we display results from a new model of Endresen et al. , where the voltage is not a dynamic variable. Here we cannot select an initial voltage independently of the initial charge, and only one limit cycle is observed.

## II Existing Models

The model of Wilders et al.  of the pacemaker activity of the rabbit sinoatrial node serves as an excellent example of the many models where the membrane potential is thought to be independent of the intracellular and extracellular charge concentrations. In that model the equation of motion for the voltage is given by (2)

$$\frac{dv}{dt}=\frac{1}{C}(i_{\mathrm{b},\mathrm{Ca}}+i_{\mathrm{b},\mathrm{K}}+i_{\mathrm{b},\mathrm{Na}}+i_{\mathrm{Ca},\mathrm{L}}+i_{\mathrm{Ca},\mathrm{T}}+i_\mathrm{f}+i_\mathrm{K}+i_{\mathrm{Na}}+i_{\mathrm{NaCa}}+i_{\mathrm{NaK}}).$$
(3)   
There are fifteen dynamic variables in that model, the voltage $`v`$, the gating variables $`d_\mathrm{L}`$, $`d_\mathrm{T}`$, $`f_\mathrm{L}`$, $`f_\mathrm{T}`$, $`x`$, $`y`$, $`h`$, $`m`$, $`p`$, and the ionic concentrations $`[\mathrm{Ca}]_\mathrm{i}`$, $`[\mathrm{Ca}]_{\mathrm{rel}}`$, $`[\mathrm{Ca}]_{\mathrm{up}}`$, $`[\mathrm{K}]_\mathrm{i}`$, $`[\mathrm{Na}]_\mathrm{i}`$. We want to determine how the long term dynamics in that model is changed when we change the initial conditions. To keep matters simple, we only change one initial condition: the initial intracellular concentration of potassium ($`[\mathrm{K}]_\mathrm{i}`$); and we study the dynamics in two dimensions only: the phase space of $`v`$ and $`[\mathrm{K}]_\mathrm{i}`$.

Figure 1 displays the two–dimensional dynamics of the model with the initial conditions (dimensions skipped)

$$\begin{array}{ccc}v=60.03\hfill & x=0.3294906\hfill & [\mathrm{Ca}]_\mathrm{i}=0.0000804\hfill \\ d_\mathrm{L}=0.0002914\hfill & y=0.1135163\hfill & [\mathrm{Ca}]_{\mathrm{rel}}=0.6093\hfill \\ d_\mathrm{T}=0.0021997\hfill & h=0.1608417\hfill & [\mathrm{Ca}]_{\mathrm{up}}=3.7916\hfill \\ f_\mathrm{L}=0.9973118\hfill & m=0.1025395\hfill & [\mathrm{K}]_\mathrm{i}=\mathrm{variable}\hfill \\ f_\mathrm{T}=0.1175934\hfill & p=0.2844889\hfill & [\mathrm{Na}]_\mathrm{i}=7.5,\hfill \end{array}$$

and $`[\mathrm{K}]_\mathrm{i}=140`$. In figure 1 (b) the point $`\mathrm{P}_1`$ denotes the $`v`$ and $`[\mathrm{K}]_\mathrm{i}`$ coordinates of the initial conditions, and the closed loop at the bottom is the limit cycle. The trajectory of the model is displayed in 1 (a) where we observe that the model spirals downwards from the point $`\mathrm{P}_1`$ to the limit cycle.

If the limit cycle in figure 1 (b) is unique it should be possible to reach it from another initial condition. Let us try an initial condition below the limit cycle, and investigate whether the model spiral up towards it. We change the initial concentration of potassium from $`140`$ to $`130`$, leaving the fourteen other initial conditions unchanged. The result is displayed in figure 2. The model does not spiral upwards to the limit cycle in figure 1, instead the model spiral downwards to a different limit cycle. We observed numerically a new limit cycle for each new initial value of $`[\mathrm{K}]_\mathrm{i}`$, implying the existence of an infinite number of limit cycles. The model’s fundamental flaw is clearly demonstrated.

## A New Model

In a new model  of the pacemaker activity of the rabbit sinoatrial node, the membrane potential is determined by (1)

$$v=\frac{FV}{C}\left\{[\mathrm{K}]_\mathrm{i}[\mathrm{K}]_\mathrm{e}+2([\mathrm{Ca}]_\mathrm{i}[\mathrm{Ca}]_\mathrm{e})+[\mathrm{Na}]_\mathrm{i}[\mathrm{Na}]_\mathrm{e}\right\},$$
(4)   
where $`q=FV\left\{[\mathrm{K}]_\mathrm{i}[\mathrm{K}]_\mathrm{e}+2([\mathrm{Ca}]_\mathrm{i}[\mathrm{Ca}]_\mathrm{e})+[\mathrm{Na}]_\mathrm{i}[\mathrm{Na}]_\mathrm{e}\right\}`$ is the charge difference, $`F`$ is Faraday’s constant and $`V`$ is cell volume. Here the ionic currents alter the concentrations which in turn alter the voltage, i.e. the physical quantities were calculated in the following order: $`\mathrm{current}i\mathrm{charge}q\mathrm{voltage}v`$. The model has five dynamic variables, the gating variables $`x`$, $`h`$ and the ionic concentrations $`[\mathrm{Ca}]_\mathrm{i}`$, $`[\mathrm{K}]_\mathrm{i}`$, $`[\mathrm{Na}]_\mathrm{i}`$, and we use the initial conditions:

$$\begin{array}{c}x=0.9165\hfill \\ h=0.0000\hfill \\ [\mathrm{K}]_\mathrm{i}=\mathrm{variable}\hfill \\ [\mathrm{Ca}]_\mathrm{i}=0.004094141\hfill \\ [\mathrm{Na}]_\mathrm{i}=18.73322695.\hfill \end{array}$$
(5)   
In this model we first notice that the initial value of $`[\mathrm{K}]_\mathrm{i}`$, due to (4), is not independent of the initial value of the voltage $`v`$. Thus changing $`[\mathrm{K}]_\mathrm{i}`$ changes $`v`$ as is always the case when charging a capacitor (1). Second we notice that a tiny change in $`[\mathrm{K}]_\mathrm{i}`$ corresponds to a large change in voltage $`v`$, since the constant $`FV/C`$ is large in most situations.

In figure 3 we have displayed the simulation results from the model of Endresen et al.  with three slightly (due to the large constant $`FV/C`$) different initial values of $`[\mathrm{K}]_\mathrm{i}`$: 130.650 (a), 130.655 (b), and 130.662 (c). In figure 3 (d) the three initial conditions $`\mathrm{P}_1`$, $`\mathrm{P}_2`$ and $`\mathrm{P}_3`$ all converge towards the same limit cycle. In an extensive numerical study we have not observed any physiological initial conditions that do not converge toward this limit cycle. In fact the same limit cycle can be reached when starting from the full equilibrium situation with equal intracellular and extracellular ionic concentrations .

## III Discussion

We have displayed numerical results from two types of mathematical models of the pacemaker activity of the rabbit sinoatrial node. The first type of model  showed an infinite number of limit cycles, the second type of model  a limit cycle that could be reached from many different initial conditions. In order to avoid the drawback with an infinite number of limit cycles seen in the first type of models, we suggest that one should not treat membrane voltage ($`v`$) as a dynamic variable. Instead one should calculate the voltage using (1), or at least select the initial conditions in agreement with (1) .

## Acknowledgments

Lars Petter Endresen was supported by a fellowship at NTNU, and has received support from The Research Council of Norway (Programme for Supercomputing) through a grant of computing time.