Quantitative Decomposition of Dynamics of Mathematical Cell Models: Method and Application to Ventricular Myocyte Models.

Mathematical cell models are effective tools to understand cellular physiological functions precisely. For detailed analysis of model dynamics in order to investigate how much each component affects cellular behaviour, mathematical approaches are essential. This article presents a numerical analysis technique, which is applicable to any complicated cell model formulated as a system of ordinary differential equations, to quantitatively evaluate contributions of respective model components to the model dynamics in the intact situation. The present technique employs a novel mathematical index for decomposed dynamics with respect to each differential variable, along with a concept named instantaneous equilibrium point, which represents the trend of a model variable at some instant. This article also illustrates applications of the method to comprehensive myocardial cell models for analysing insights into the mechanisms of action potential generation and calcium transient. The analysis results exhibit quantitative contributions of individual channel gating mechanisms and ion exchanger activities to membrane repolarization and of calcium fluxes and buffers to raising and descending of the cytosolic calcium level. These analyses quantitatively explicate principle of the model, which leads to a better understanding of cellular dynamics.

Introduction

Mathematical modelling has been an effective method in physiology for precise and comprehensive understanding of the dynamic behaviour of cells. A number of mathematical cell models have been developed, and recent models of cardiac cells [1-5] have been more detailed and thereby complicated by including multiple cellular functions to explain new experimental findings. Conventionally, these models have been used to simulate wet experiments. In contrast with wet experiments, an accurate and more complete set of experimental data can be obtained by numerical simulation. Additionally, mathematical models enable simulation experiments that are otherwise impracticable, such as a pure and complete blockade of an ion channel or a perfect control of the intracellular composition. Despite the success of simulation, such conventional simulation is insufficient to achieve the full potential of mathematical cell models. Since the whole mechanisms of each model dynamics are explicitly defined in mathematical expressions, models potentially enable quantitative clarification of their detailed behaviour, which leads to a better understanding of cellular dynamics.

Each of mathematical cell models is generally formulated as a system of ordinary differential equations (ODE) with respect to time. The ODE model variables interact with each other either directly or indirectly and vary simultaneously. In order to elucidate the causes and results of this interaction, inspection of model equations is essential but difficult for detailed models due to complicated interdependences of variables. To overcome this difficulty, mathematical approaches are required. One such approach applicable to mathematical cell models is bifurcation analysis, which is used to investigate qualitative changes in a system of equations by smooth changes in parameter values. More specifically, the bifurcation analysis can determine whether a model converges, diverges, or oscillates depending on the parameter values. For instance, Kurata and his collaborators [6-12] have applied the bifurcation analysis to mathematical models for understanding the oscillatory phenomena in ventricular and sinoatrial node cells. The singular perturbation method of asymptotic analysis is a method for inspection of the dynamic behaviour of mathematical models. In this method, variables are divided into fast and slow ones, and steady states of a model in regarding the slow variables as parameters are traced in time. Analysis based on this method can explain dynamic change in characteristics, e.g. membrane excitability of cardiac cells [13-16]. These methods can answer why a model has its behaviour.

Another fundamental question in model dynamics is how much each model component affects the model behaviour. In physiological experiments, the most conventional approach for examining contribution of a cellular component is activation or inhibition of a target function using agonists, blockers or knockout of the corresponding gene. The same kinds of methods have been also applied to many simulation studies by altering the corresponding parameter values. However, the interpretation of results of these methods for estimating contribution of a component in physiological condition is extremely difficult in most cases. Since a modification to a component secondarily causes changes in other components which also affect the target function, the resultant change in the function cannot be considered as a sole effect of the modulated component but a mixed effect of the other components. To overcome this difficulty, Clewley et al. [17, 18] have developed ‘dominant scale method’, and Cha et al. [19] ‘lead potential analysis’. However, their methods are limited to analyses of cellular membrane potential.

In this study, a numerical method is introduced for quantitatively decomposing dynamics of mathematical cell models. This method is applicable to analysis of every model variables, and able to evaluate contributions of individual model components to the dynamics of a variable. Firstly, the mathematical definition of the proposed method is presented in this article. Then, applications of the method to action potential and calcium transient of ventricular myocyte models are illustrated.

Method

This section provides the mathematical definition of a novel index for decomposed dynamics of an object variable in an ODE model, following to introduction of a concept ‘instantaneous equilibrium point.’ For a time-dependent variable v, v ( denotes the value of v at t, and denotes the first derivative of v with respect to t, i.e., dv/dt. A variable v is called a differential variable if is explicitly defined in a model.

Consider the dynamics of an object differential variable x. Let the derivative of x where y is a vector of all the other differential variables.

By introducing an auxiliary time-dependent function such that at t for any value of x, for simplicity, Eq 1 at a certain time t can be written as By linearizing Eq 2 with the first-order Taylor series of at x (, where stands for the partial derivative of f with respect to x at (x (, y (). Eq 3 expresses the tangent line to at x (, or to f(x, y) at the current value of x in assuming that y is constant with its value at t . The x-intercept of the tangent line is . Consequently, the linearized Eq 1 at t with respect to x is obtained: Thus, is the fixed point of Eq 5, which is stable if is negative or unstable if positive. Namely, x is naively attracted to or repelled from . Hence, is referred to as the instantaneous equilibrium point of x, which represents the trend in x at some instant. Since is time-dependent, the orbit of x in Eq 1 is additionally affected by temporal change in . Therefore, temporal change in can be regarded as representing active dynamics of x.

To decompose the active dynamics of x, consider an effect of temporal change in a variable y of y on a variation in at a certain moment. A sensitivity k of the rate of change in to the rate of change in y is mathematically expressed as where f stands for the second-order mixed derivative of f with respect to x and y, and f for the partial derivative of f with respect to y. Here, y dynamic of x is defined by the sensitivity k weighted by , the current rate of change in y at the moment: Similarly, x dynamic of x itself is defined: where f stands for the second-order partial derivative of f with respect to x.

For the chain rule, Consequently, c is component of the rate of change in , that is, the decomposed dynamics of x. The sign of c , which is determined by the both signs of k and , indicates whether a temporal change in y induces x to increase in the positive case or decrease in the negative case. Note that c represents an effect of the nonlinearity of .

Application

Models and Methods

The present method is demonstrated using the two mathematical ventricular myocyte models. One is a guinea pig model proposed by Takeuchi et al. [1], which represents membrane excitation, ion homeostasis, excitation-contraction coupling, volume regulation, and the balance of ATP production and consumption. This model formulates the spermine block and the magnesium block of I K1, which is a major outward current of ventricular myocytes. In the original model, while the open probability y of the spermine block is expressed as a differential equation, the open probability f of the magnesium block is defined as a steady-state expression. For that expression a differential equation is substituted in this application (see S1 Appendix), in order to evaluate the contributions of I K1 through both the magnesium and spermine blocks. As a result, the number of differential variables of the modified model is 51. The other is a human ventricular myocyte model published by Priebe and Beuckelmann [5], which includes action potential generation and calcium dynamics. The number of its differential variables is 22.

A simulation program is coded in C to invoke an ODE solver, CVODE (Lawrence Livermore National Laboratory, Livermore, CA), in which backward differentiation formulae are used and the time step is adaptively controlled. The first-order partial derivatives of the derivative functions of differential variables are obtained with automatic formula differentiation. The second-order partial derivatives are calculated with numerical differentiation. The stimulation protocol applied to each model follows the corresponding original paper [1, 5]. For the Takeuchi model, the stimulation current is injected for 2 ms at 50 ms after the start of the simulation. For the Priebe model, the stimulation is applied for 3 ms at 0 ms. The values of parameters and initial values of variables used in the simulation are identical to the original models.

In order to analyse the mechanisms of action potential generation and calcium transient, which are central functions in ventricular myocyte, object variables in the present analysis are the membrane potential and the intracellular calcium concentration. For the Takeuchi model, the amount of cytosolic calcium n(Ca)i replaces the intracellular calcium concentration, because the concentration is not a differential variable but defined as n(Ca)i divided by the cellular volume. For the Takeuchi model, simulated time courses of the membrane potential V and n(Ca)i to be analysed are shown with black lines in Figs 1 and 2, respectively. The Takeuchi model reproduces a typical ventricular action potential with a long plateau, which is formed by inward and outward ion currents through cellular membrane. The major inward currents are I Na, I CaL, I CaT and I NaCa, and the outward currents are I K1 and I l(Ca). An intracellular calcium transient is determined in the Takeuchi model by transmembrane Ca2+ currents (I CaL, I CaT and I NaCa) and Ca2+ fluxes across the sarcoplasmic reticulum (SR) membrane (I RyR and I SERCA). These ion currents and fluxes are shown with black areas in Fig 3. Differential variables for open probabilities of channel gates, state variables of transporters, and calcium concentration in SR are also plotted with solid lines. In a similar manner, for the Priebe model, simulated time courses of the membrane potential V and the intracellular calcium concentration [Ca2+]i are plotted with black lines in Figs 4 and 5, respectively. All the time courses of model differential variables are identical to those in [1] and [5].

Fig 1

Time course of membrane potential (black) and its instantaneous equilibrium point (red) of Takeuchi model.

Fig 2

Time course of amount of cytosolic calcium (black) and its instantaneous equilibrium point (red) of Takeuchi model.

Fig 3

Major currents and differential variables simulated using Takeuchi model.

(A) I Na, fast sodium current; I Na p(AP), open probability of I Na voltage gate. (B) I K1, inward rectifier potassium current; I K1 f O, open probability of I K1 magnesium gate; I K1 y, open probability of I K1 polyamine gate. (C) I CaL, L-type calcium current; I CaL p(AP), open probability of I CaL voltage-dependent gate; I CaL p(U), open probability of I CaL calcium gate; I CaL y, open probability of I CaL ultra-slow gate. (D) I CaT, T-type calcium current; I CaT y 1, open probability of I CaT activation gate. (E) I NaCa, Na+/Ca2+ exchange current. (F) I , Ca2+-activated background cation current. (G) I RyR, ryanodine receptor channel current; I RyR p(Open), open probability of I RyR gate. (H) I SERCA, SR Ca2+ pump current; I SERCA y, probability of the conformation state with the Ca2+-binding sites onto the SR side. (I) [Ca2+]rel, calcium concentration in SR release site.

Fig 4

Time course of membrane potential (black) and its stable and unstable instantaneous equilibrium point (red and blue, respectively) of Priebe model.

Fig 5

Time course of intracellular calcium concentration (black) and its instantaneous equilibrium point (red) of Priebe model.

Dynamics of Action Potential

For the Takeuchi model, a red line in Fig 1 represents the instantaneous equilibrium point (IEP) of V at each time defined in Eq 4, where V , and the partial derivative of function with respect to V are substituted for x, and f , respectively. Specifically, the IEP of V , is expressed as The IEP is higher than V during depolarization phase but lower during repolarization phase, that is, V is always attracted to its IEP, because the IEP is stable or f in Eq 5 is negative throughout the simulation period. The orbit of the IEP intersects with the time course of V exactly at the peak of V (∼ 53 ms) for Eq 5. Note that IEP during 50–52 ms is not effective to analysis due to the external current injection. Along with V time course and its IEP on the same time base in the upper panel, each variable dynamic of V with a significant value is plotted in the lower panels of Fig 6, separated into four sequential phases of the cardiac action potential; the last 0.5 ms of rising phase of the action potential after termination of the triggering current injection (phase 0, Fig 6A), the initial repolarizing phase of 4 ms (phase 1, Fig 6B), the major plateau phase (phase 2, Fig 6C), and the final repolarization (phase 3, Fig 6D). Note that, for a differential variable v, v dynamic of V is expressed as corresponding to Eq 7. Similarly for the Priebe model, a red line in Fig 4 represents the stable IEP of V, and a blue line during 323–365 ms represents the unstable IEP of V. Each variable dynamic of V with a significant value is plotted in the lower panels of Fig 7, separated into four sequential phases; the period including the V peak of 2 ms after termination of the triggering current injection (phase 0, Fig 7A), the initial repolarizing phase (phase 1, Fig 7B), the major plateau phase (phase 2, Fig 7C), and the final repolarization (phase 3, Fig 7D). IEP during 0–3 ms is not analysable due to the stimulation current.

Fig 6

Dynamics of the membrane potential V of Takeuchi Model after stimulus offset at different phases (A–D) with different scales, together with V (black) and the instantaneous equilibrium point (red) at the top of each figure on the same time base.

I NaCa denotes the sum of I NaCa p(E1total), p(I1) and p(I2) dynamics of V . A notch of V dynamic around 185 ms and a tiny fluctuation of I Na p(AP) dynamic around 235 ms in D are errors associated with division by nearly zero in numerical differentiation.

Fig 7

Dynamics of the membrane potential V of Priebe model after stimulus offset at different phases (A–D) with different scales, together with V (black) and the instantaneous equilibrium point (red) at the top of each figure on the same time base.

An increase in the IEP of the membrane potential prompts an increase in the membrane potential, and vice versa. Therefore, a variable v induces a membrane depolarization if the v dynamic of V (or V) is positive, or a repolarization if negative. The sign of the v dynamic of V is determined by both k = ∂V /∂v and the derivative of v. Since an inward current itself induces an increase in V , k for an open probability of an inward current is positive. Therefore, a decrease in the open probability, or its negative derivative, induces a decrease in the IEP of V and contributes to repolarization, while an increase in the probability contributes to depolarization. For an outward current, the effect is opposite.

For the Takeuchi model, Fig 6A shows that the fast depolarization in phase 0 is exclusively attributed to the V dynamics of V itself. The details are explained later. In phase 1 shown in Fig 6B, I Na p(AP), the open fraction of I Na voltage gates is the dominant negative dynamic of V . On the other hand, as shown in Fig 7A for the Priebe model, the IEP of V is already descending at the termination time of the stimulation current, and V simply increases towards the IEP. After the intersection of V and its IEP (∼ 3.5 ms), V decreases in pursuit of the IEP until 10 ms (Fig 7B). The main contributor to the decrease in the IEP is m, the activation gate of I Na (Fig 7A). These results for both models agree with an accepted view that inactivation of I Na voltage gate mainly causes the initial negative shift in V . The reason for the large negative dynamic is that the open fraction is decreasing due to rapid inactivation while the reversal potential of I Na is more positive than the IEP. For the Priebe model, V increases again after its initial negative shift, because the IEP exceeds V again (Fig 7B). The major contributors to the increase of the IEP are r and t, the activation gate and the inactivation gate of I to, respectively. This also agrees with a physiological view.

In phase 2 shown in Fig 6C for the Takeuchi model, the major contributors to membrane repolarization are I CaL p(U) (open probability of I CaL calcium gates), and n(Ca)i subsequently replaced by I NaCa (Na+/Ca2+ exchange current). I CaL p(U) dynamic of V is negative for all phase 2 and least during 65–110 ms. Decreasing of I CaL p(U) (Fig 3C) by gradual Ca2+-dependent inactivation results in a reduction of I CaL, which provides a significant inward current. This analysis shows explicitly that an inward current, which has been generally thought to induce depolarization, contributes to a repolarization when its open probability is decreasing. Two other variables in I CaL, voltage-dependent inactivation (I CaL p(AP)) and ultra-slow inactivation (I CaL y), also have negative dynamic values. Therefore, closing of I CaL gates induces the repolarization in phase 2. From the opposite side, this is consistent with a physiological view that I CaL has major contribution to forming the plateau phase of action potential. The same is true for the Priebe model. As shown in Fig 7C, f, the inactivation gate of I Ca has negative and least dynamic of V until 175 ms. For the Takeuchi model, the n(Ca)i dynamic of V is least during the early period of the phase 2 about 58–66 ms, but becomes positive after the peak of n(Ca)i shown in Fig 2. These indicate that the rise in n(Ca)i has a large repolarizing effect on the Takeuchi model. It is mainly attributable to an activation of I l(Ca) (Fig 3F). Since I l(Ca) is defined as directly depending on V and n(Ca)i without any state variables in the model, its effect on the dynamics of V is not explicitly shown, and included in n(Ca)i dynamics of V in this method. In Fig 7C for the Priebe model, X and X , the activation gate of I Ks and I Kr respectively, have relatively small negative dynamics of V. In Fig 6C for the Takeuchi model, dynamics of V relating to I Ks or I Kr are not shown, because they are negligibly small. These results indicate that I Ks and I Kr make effects on the level of membrane potential in phase 2 but rather small dynamic effects to induce repolarization.

In the classical view, final repolarization towards the resting potential in phase 3 is mainly attributable to increase in I K1, decrease in I CaL, and deactivation of I NaCa. The results in Fig 6D show general agreement with the view, and give additional insights into the Takeuchi model. Voltage-dependent removal of Mg2+ block of I K1 (I K1 f O) induces the repolarization, whereas closing of voltage-dependent gates (I K1 y) prevent it slightly at late phase 3. Similarly, deactivation of I CaL voltage gates (I CaL p(AP)) induces the repolarization and opening of Ca2+-dependent gates (I CaL p(U)) hinders it. Finally, the negative I NaCa dynamic of V indicates a repolarizing effect of I NaCa, contrary to the fact that its inward current is increasing until ∼ 225 ms (Fig 3E). The negative I NaCa dynamic, which denotes the sum of I NaCa p(E1total), p(I1) and p(I2) dynamics, results from a rearrangement of molecular states of the exchanger to decelerate the exchange of intracellular Ca2+ for extracellular Na+. On the other hand, an inward current of I NaCa is determined not only by the molecular states but also by driving force of Ca2+ influx through the exchanger, which is enlarged by membrane repolarization. Thus, I NaCa affects the repolarization and is affected by it at the same time. The present method can discriminate the former active effect from the latter passive effect (see Discussion).

As shown in Fig 6A for the Takeuchi model, the V dynamic of V itself is dominant in phase 0. It indicates that the increase in V induces further depolarization. In this period, massive I Na, which is formulated with Goldman-Hodgkin-Katz flux equation, gives high nonlinearity to the derivative function of V . For an instance, the derivative of V at 52.05 ms is shown in Fig 8. Due to the nonlinearity, a positive change in V simply shifts the x-intercept, or the IEP of V , in the positive direction. In Fig 7D for the Priebe model, the V dynamic of V itself is dominant after 310 ms of the final repolarization phase. Fig 9 shows the derivatives of V at 257 ms and 340 ms. When V is greater than the local maximum of the derivative (Fig 9A), V is headed for its stable IEP, and a negative perturbation to V makes the IEP less because of the nonlinearity of the derivative. The IEP approaches to infinity as V gets close to the local maximum. When V exists between the local maximum and minimum (Fig 9B), V tends to move in the opposite direction of its unstable IEP, and a negative change in V shifts its IEP in the negative direction. Therefore, a reduction in V decreases the IEP, and this induces repolarization automatically in this phase. Note that enormous negative and positive dynamics of V balances around switching points of the stability of the IEP of V, because the derivative of the IEP is extremely large (Eq 9).

Fig 8

Derivative function of V of Takeuchi model at 52.05 ms (thick curve) assuming that all the other differential variables are constant.

× is the instantaneous equilibrium point of V , which is the x-intercept of the tangent line (thin line) at the value of V at that time (+).

Fig 9

Derivative function of V of Priebe model at 257 ms and 340 ms (thick curve in A and B, respectively) assuming that all the other differential variables are constant.

× is the instantaneous equilibrium point of V, which is the x-intercept of the tangent line (thin line) at the value of V at that time (+).

Dynamics of Calcium Transient

Another practice of the present method is on calcium transient, which is a transient elevation of intracellular Ca2+ concentration. For the Takeuchi model, the amount of cytosolic calcium n(Ca)i and its IEP are plotted with a black and red lines in Fig 2, respectively. The IEP leads the time course of n(Ca)i because the IEP is stable throughout the simulation period. In the Takeuchi model, the calcium transient is determined by balance among three mechanisms; (1) the net Ca2+ flux across the surface membrane via L-type and T-type Ca2+ channels (I CaL, I CaT), Na+/Ca2+ exchangers (I NaCa) and plasma membrane calcium pumps (I PMCA), (2) the net Ca2+ flux across the SR membrane via ryanodine receptors (I RyR) and calcium pumps (I SERCA), and (3) Ca2+-binding to cytosolic Ca2+ buffer proteins. For measuring contributions of these factors, each variable dynamic of n(Ca)i with a significant value are plotted along with total dynamics, the sum of all the dynamics, in two periods; rising and plateau phases of n(Ca)i in Fig 10A, and falling phase in Fig 10B. The total dynamics of n(Ca)i is equivalent to the derivative of the IEP of n(Ca)i for Eq 9. For the Priebe model also, intracellular calcium concentration [Ca2+]i and its stable IEP are plotted with a black and red lines in Fig 5, respectively. Individual dynamics of [Ca2+]i and the total dynamics are shown in two phases; upward phase of [Ca2+]i in Fig 11A, and downward phase in Fig 11B. Similarly to the analysis of the membrane potential, a positive dynamic of n(Ca)i or [Ca2+]i indicates an effect to vary the IEP in the positive direction, and vice versa.

Fig 10

Each dynamic and total dynamics of amount of cytosolic calcium, n(Ca)i of Takeuchi model after stimulus onset at two phases (A, B) with different scales, together with n(Ca)i (black) and the instantaneous equilibrium point (red) at the top of each figure.

I NaCa stands for the sum of I NaCa p(E1total), p(I1) and p(I2) dynamics of n(Ca)i, and troponin stands for the sum of [T], [TCa], and [TCa*] dynamics of n(Ca)i.

Fig 11

Each dynamic and total dynamics of intracellular calcium concentration, [Ca2+]i of Priebe model after stimulus onset at two phases (A, B) with different scales, together with [Ca2+]i (black) and the instantaneous equilibrium point (red) at the top of each figure.

In the general physiological view, the rising phase of the cytosolic calcium level is attributable to calcium-induced calcium release (CICR), in which a slight increase in intracellular Ca2+ by the activation of I CaL and I CaT further triggers a massive Ca2+ release from sarcoplasmic reticula into cytosol through ryanodine receptors (RyRs). In Fig 10A for the Takeuchi model, the total dynamics of n(Ca)i has the first very brief positive peak during rapid depolarization of the cellular membrane. This peak is formed by I CaL p(AP) dynamic and I CaT y 1 dynamic of n(Ca)i and additionally by the first sharp peak of I RyR p(Open) dynamic, which is caused by an activation of RyRs proportional to Ca2+ flux through I CaL. These positive dynamics drive n(Ca)i in the positive direction slightly. The second wave of the total dynamics is mainly attributed to I RyR p(Open) via an increase in cytosolic Ca2+ itself caused by the first peak of the total dynamics. These analysis results of the Takeuchi model are in good agreement with the view of the CICR mechanism. However, subsequent results reveal an important difference. I RyR p(Open) begins to decrease at 57.8 ms (Fig 3G) during the rise in the IEP of n(Ca)i due to time-dependent closure of RyRs. This decrease makes I RyR p(Open) dynamic negative and hinders the rise in n(Ca)i. In addition, the Ca2+ release via RyRs rapidly decreases [Ca2+]rel, total calcium concentration in SR release site (Fig 3I). This decrease also hampers the increase in n(Ca)i via reduction in the Ca2+ concentration gradient across the SR membrane to drive I RyR, as indicated by negative [Ca2+]rel dynamic in Fig 10A. Alternatively, the prominent contributor in the Tekeuchi model to the increase in the IEP of n(Ca)i after 55 ms is I SERCA, which is to reduce n(Ca)i by pumping cytosolic Ca2+ into SR though. In a rough explanation, deactivation of I SERCA decreases Ca2+ flux from cytosol to SR and so induce to increase n(Ca)i. To be precise, I SERCA y is the fraction of carrier proteins exposing their Ca2+-binding sites toward SR lumens and negatively represents a capability of I SERCA to transfer calcium ions from cytosol. I SERCA y dynamic is positive after 52 ms until 77 ms because of the increase in I SERCA y, although I SERCA is increasing until 56 ms and decreasing still after 77ms (Fig 3H) due to passive dependency of I SERCA on cytosolic Ca2+ concentration. Additionally, the results show that troponin, which is a cytoplasmic calcium buffer protein, partially contributes to the rise in the IEP of n(Ca)i until 104 ms (Fig 10A and 10B) by a decrease in concentration of calcium-free troponin, or Ca2+ binding sites.

The negative total dynamics of n(Ca)i shown in Fig 10B facilitates the last calcium decrease. I RyR p(Open) dynamic is continuously negative since the late rising phase until 162 ms. Additionally, negative I SERCA y dynamic via a decline in I SERCA y and negative troponin dynamic by the opposite of the aforementioned effect account for a large portion of the total dynamics of n(Ca)i. In contrast, I NaCa dynamic of n(Ca)i is positive due to increased proportion of inactive conformation states of the exchanger, whereas an inward current of I NaCa (Fig 3E) to carry Ca2+ out of cell is increasing until the end of action potential (∼ 225 ms) because of an increase in the voltage-dependent transport rate of I NaCa. An effect of the voltage-dependency on n(Ca)i is represented by a major part of negative V dynamic. Therefore, overall contribution of I NaCa to dynamic changes of cytosolic calcium level is relatively small.

For the Priebe model, the mechanism of calcium transient is rather simple, because the calcium release and uptake of SR and I NaCa are passively formulated without differential variables. As shown in Fig 11, the elevation of [Ca2+]i triggerd by membrane depolarization is attributed to d, the activation gate of I Ca and partially V, as is the case in the Takeuchi model. Although the IEP of [Ca2+]i descends after 12 ms caused by the deactivation of I Ca (f), the IEP is sufficiently high to raise [Ca2+]i. The time course of the IEP has a plateau around 200 ms, because d and f dynamics of [Ca2+]i come to zero and balance with each other. Finally, the main contributors to a decrease in [Ca2+]i are closing of the activation gate of I Ca (d) and voltage dependency of the transmembrane calcium flux (V).

Summary of Analyses

Analysis results in Figs 6 and 7 confirm that the Takeuchi model and the Priebe model provide generally good agreement with the classic qualitative view on the ionic mechanisms of the action potential generation in ventricular myocytes. On the other hand, these analyses exhibit that inward currents can contribute to the repolarization via their decrease, in contrast to a classic physiological concept. In the plateau phase, inactivation of Ca2+ channels prompts the repolarization. Additionally, these analyses provide quantitative insights into the contribution of each gating mechanism and ion exchanger activity. For instance, I Ks and I Kr have only small dynamic effects to induce membrane repolarization. In the last phase of membrane repolarization, Mg2+ block and voltage-dependent gates of I K1 contradictorily affect the repolarization. These insights into the dynamic mechanisms of the action potential generation will be useful for controlling the action potential.

Analyses on calcium transient of the Takeuchi model show good agreement with the general physiological view of the CICR mechanism in the early phase of [Ca2+]i elevation, but also provide an important difference. [Ca2+]i rising in the late part is attributable to the reduced capability of SERCA to take calcium ions into SR. The analyses results also show that the decay phase of calcium transient on the Takeuchi model is attributable to decrease in calcium release from SR through I RyR initially and to I SERCA and troponin thereafter.

Discussion

Feature of the Proposed Method

This article presents a novel analytical method, which is applicable to a complicated mathematical cell model with many variables, for decomposing and quantifying contributions of individual variables to dynamic change in a variable of interest at each moment during simulations. Toward this end, the method employs a representative point of the instantaneous trend in the object variable, ‘instantaneous equilibrium point’, in which the temporal change characterises dynamics of the object variable. The contributions of variables are evaluated with a quantitative index, ‘dynamic of the object variable’, of which value is a variable component of a variation in the instantaneous equilibrium point. Applications of the present method to a ventricular myocyte model demonstrate its capability to mathematically quantify dynamic interactions among variables of a mathematical model and thereby explicate principle of the model.

Comparison with Other Methods

A conventional method applied to mathematical cell models is sensitivity analysis. In a typical application manner of sensitivity analysis, variations of simulated results are observed against variations of model parameters. This analysis is useful for examining comprehensive influence of chemical or genetic modifications to cellular components. However, estimating contributions of cellular components to cellular dynamic behaviour is difficult by using this analysis, because a modification to a model parameter causes cascade effects on cellular dynamics. For an instance, a decrease in a parameter for the amplitude of Ca2+ current results in not only a decrease in the Ca2+ current but also a certain decrease in calcium level, and the decrease in calcium level affects Ca2+-dependent processes, which make further effects. In contrast, the present method is able to quantify contributions of cellular components to cellular dynamics in the intact situation. This is a distinguishing feature of the present method.

A method developed to analyze contribution of current components to an action potential is ‘lead potential analysis’ proposed by Cha et al [19]. The present method can be regarded as a generalized version of the lead potential analysis. The lead potential (V ) in the Cha’s method [19, 20] is almost equivalent to the instantaneous equilibrium point of the membrane potential, and contribution (c ) or relative contribution (r ) corresponds to the dynamic of the membrane potential. The results in Fig 6 of this article are similar to their analysis results of the Takeuchi model in Fig. 5 in [19]. Compared to the lead potential analysis, the present method has improvements and advantages in several aspects. Firstly, the c or r in Cha’s method is procedurally defined, in contrast to the mathematically well-defined present method. Next, the lead potential analysis has a limitation on treatment of ion transporters as discussed in [19]. To calculate the lead potential, the membrane current system through ion channels and transporters are approximately described with an equivalent electrical circuit. While each channel current is straightforwardly expressed in the circuit as a pair of a battery for the corresponding equilibrium potential and a variable conductance to express its dynamic gating, ion transporters are approximated as current sources [19, 21, 22] or ohmic channels [20] depending on the studies. As a result of the approximation, the contribution of a transporter includes not only effect of its activity on the membrane potential but also passive effect of the membrane potential on the transporter. Since the present method mathematically properly processes all components in the same manner, a dynamic of the membrane potential only expresses active effect. For instance, the dynamic effect of conformational changes of I NaCa on membrane potential is successfully discriminated in I NaCa dynamic of V in Fig 6, in contrast to r of I NaCa in Fig. 5 in [19]. Finally, in the lead potential analysis, selection of current components of which r are to be calculated is arbitrary and demanding, because the components should be selected so that the summation of r of all the selected components is exactly 1. In the present method, a set of variables whose dynamic are calculated is automatically and uniquely determined on the set of all differential variables. For instance, I l(Ca) dynamic is not defined in the application of the present method to the Takeuchi model because I l(Ca) is formulated with no state variable, whereas results in [19] include the contribution of I l(Ca) to the membrane potential.

Another numerical method for analysis of action potential dynamics is ‘dominant scale method’ proposed by Clewley et al [17]. with the main aim of model reduction. This method introduces the instantaneous asymptotic target, x ∞ such that f(x ∞, y () = 0 at t for , or a fixed point of Eq 2 in the present article. Then, ∂x ∞/∂y is utilized as an index of dominance, and ∣(∂V ∞/∂s)(ds/dt)∣ is referred to as the instantaneous rate of influence [18], where V is the membrane potential and s is a gating variable. The dominant scale method is identical to the application of the present method to the membrane potential if dV/dt is linear with regard to V, as in the case of Hodgkin-Huxley model shown in [18]. The primary difference of the present method from the dominant scale method is linearization of the derivative function around the current values of the object variable. Whereas the instantaneous asymptotic target is a global fixed point of a system, the instantaneous equilibrium point is to represent the instantaneous trend in an object variable as described in Method section. This is an advantage of the present method. Fig 9A displays a situation where the instantaneous asymptotic target, which is the intersection of the derivative curve and the x-axis (∼ −75 mV), is far from the instantaneous equilibrium point (∼ −44 mV). Although the derivative is zero at the instantaneous asymptotic target, the direction that V is heading in at that moment is the instantaneous equilibrium point. Another advantage of the present method is the uniqueness of the instantaneous equilibrium point. There always exists a unique instantaneous equilibrium point of an object variable of a system except if f is exactly equal to zero, whereas there can exist no or multiple instantaneous asymptotic targets. In fact, the derivative of V of the Priebe model has three instantaneous asymptotic targets in the plateau phase.

Limitations and Applications

In the present method, objects of analysis are limited to differential variables. This property is both an advantage and a disadvantage. Since the dynamics of an object variable are not defined with respect to variables of which values are passively determined by other variables, passive components are automatically excluded in analysis. However, in some models probabilities of gate opening and conformation states are not formulated as differential equations but as regular equations assuming instantaneous equilibria for a reason such as computational simplicity or numerical stability. In order to consider the corresponding dynamic effects in such cases, model equations require to be transformed into differential equations in the same manner as I K1 f of the Takeuchi model.

The present method is applicable to recent and complicated ventricular myocyte models [23, 24], although rather simple models analysed in this paper are suitable for the initial demonstration of the proposed method. Moreover, many computational models have been published for various physiological phenomena such as electrical activity, signal transduction, and metabolism. Applications of the present method are not limited to cardiac cell models but possible to any kind of models defined as a system of ordinary differential equations. On the other hand, the present method is not applicable to models that utilize formalisms other than ODE, such as cardiac cell models including detailed spatial-temporal expressions of calcium cycling.

The present method will be useful for analyses of transient phenomena and oscillatory phenomena such as afterdepolarizations of cardiac myocytes, because individual dynamics can be evaluated at any instant of a process. Control of arrhythmia is a potential application of the present analysis approach. Furthermore, the present method can analyse dynamics of a variable that contributes to dynamics of another variable successively. To take the application to the Takeuchi model as an instance, since n(Ca)i has a primary effect on decrease in V during the rising phase of n(Ca)i and I SERCA y is the main contributor to the rise in n(Ca)i, it is quantitatively clarified that the activity of I SERCA indirectly affects the membrane repolarization. Such consecutive analyses will help quantitative explanations of pathways of physiological phenomena.

Modified IK1 Magnesium Blockade Model.

(PDF)

Click here for additional data file.