Fast-slow asymptotics for a Markov chain model of fast sodium current

TL;DR

This paper investigates the use of fast-slow asymptotic methods to reduce computational stiffness in Markov chain models of cardiac sodium channels, aiming to improve simulation efficiency.

Contribution

It introduces a novel application of asymptotic analysis to simplify Markov chain models of ionic channels in cardiac electrophysiology.

Findings

Demonstrates feasibility of asymptotic reduction

Reduces computational stiffness in models

Provides insights into sodium channel dynamics

Abstract

We explore the feasibility of using fast-slow asymptotic to eliminate the computational stiffness of the discrete-state, continuous-time deterministic Markov chain models of ionic channels underlying cardiac excitability. We focus on a Markov chain model of the fast sodium current, and investigate its asymptotic behaviour with respect to small parameters identified in different ways.

Tables

Table 1. Table 1: Ranges of indices used in the text, unless explicitly stated otherwise.

index	values	corresponds to
$i,j,j_1,j_2$	$\iota ,\dots ,n$ 111 Here $\iota =1$ for Section II and $\iota =0$ from the next section on.	all eigenvalues of Jacobian (J.)
$k,p$	$\iota ,\dots ,m$ 111 Here $\iota =1$ for Section II and $\iota =0$ from the next section on.	zero eigenvalues of J.
	$0$ 222 This is used starting from Section III.	autonomous time direction
$q,q^{\prime }$	$1,\dots ,n$ 222 This is used starting from Section III.	all eigenvalues of Markov chain (M.C.)
$r,r^{\prime }$	$1,\dots ,m$ 222 This is used starting from Section III.	zero eigenvalues of M.C. ${\lambda }_r=0$
$\ell$	$m+1,\dots ,n$	non-zero eigenvalues of J. or M.C.

Equations

$$\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}=\mathbf{f}(\mathbf{u})+\epsilon \mathbf{h}(\mathbf{u})$$

$$\mathbf{u}=\mathbf{U}(\mathbf{a})+\epsilon \mathbf{v}(t)$$

$$\mathbf{v}(t)=\sum _{\ell }{b}_{\ell }{\mathbf{V}}_{\ell }(\mathbf{a}),$$

$$\hat{F}{\mathbf{V}}_i={\Lambda }_i{\mathbf{V}}_i,$$

$${\mathbf{V}}_k=\frac{\partial \mathbf{U}}{\partial a_k}.$$

$$\frac{1}{\epsilon }\frac{\mathrm{d}{a}_k}{\mathrm{d}t}=$$

$$\frac{\mathrm{d}{b}_{\ell }}{\mathrm{d}t}=$$

$${\mathcal{F}}_i(\mathbf{a},\mathbf{b})=$$

$$\text{mbox}+\frac{1}{\epsilon }\sum _k\frac{\partial {\mathbf{W}}_i^T}{\partial a_k}\frac{\mathrm{d}{a}_k}{\mathrm{d}t}\mathbf{v},$$

$$\hat{H}(\mathbf{U})={\frac{\partial \mathbf{h}}{\partial \mathbf{u}}|}_{\mathbf{u}=\mathbf{U}}.$$

$$b_{\ell }=-\frac{{\mathbf{W}}_{\ell }^T\mathbf{h}(\mathbf{U})}{{\Lambda }_{\ell }}+\mathcal{O}(\epsilon ),$$

$$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\hat{A}(t)\mathbf{x}.$$

$$\mathbf{u}=\left[\begin{array}{c}\sigma \\ \mathbf{x}\end{array}\right]$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}\sigma \\ \mathbf{x}\end{array}\right]=\left[\begin{array}{c}1\\ \hat{A}(\sigma )\mathbf{x}\end{array}\right].$$

$$\hat{A}=\frac{1}{\epsilon }{\hat{A}}_f+{\hat{A}}_s.$$

$${\hat{A}}_f(\sigma ){\mathbf{\kappa }}_q(\sigma )=$$

$${\hat{A}}_f^T{\mathbf{\rho }}_q={\lambda }_q{\mathbf{\rho }}_q,{{\mathbf{\rho }}_q}^T{\mathbf{\kappa }}_{q^{\prime }}={\delta }_{q,q^{\prime }}.$$

$$\frac{\mathrm{d}{{\mathbf{\rho }}_i}^T}{\mathrm{d}\sigma }{\mathbf{\kappa }}_j=-{{\mathbf{\rho }}_i}^T\frac{\mathrm{d}{\mathbf{\kappa }}_j}{\mathrm{d}\sigma }.$$

$$\frac{\mathrm{d}\sigma }{\mathrm{d}\tau }=$$

$$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}\tau }=$$

$$\mathbf{u}=\left[\begin{array}{c}\sigma \\ \mathbf{x}\end{array}\right],\mathbf{f}=\left[\begin{array}{c}0\\ {\hat{A}}_f(\sigma )\mathbf{x}\end{array}\right],\mathbf{h}=\left[\begin{array}{c}1\\ {\hat{A}}_s(\sigma )\mathbf{x}\end{array}\right].$$

$$\mathbf{U}(\sigma ,\mathbf{a})=\left[\begin{array}{c}{a}_0\\ {\mathbf{U}}^x(\mathbf{a})\end{array}\right]$$

$${\mathbf{U}}^x(\mathbf{a})=\sum _r{a}_r{\mathbf{\kappa }}_r.$$

$$\hat{F}(\mathbf{U})=\left[\begin{array}{cc}0& 0{}^T\\ \frac{\mathrm{d}{\hat{A}}_f}{\mathrm{d}\sigma }{\mathbf{U}}^x& {\hat{A}}_f\end{array}\right].$$

$${\mathbf{V}}_i=\left[\begin{array}{c}{\mu }_i\\ {\mathbf{V}}_i^x\end{array}\right].$$

$${\Lambda }_i{\mu }_i$$

$${\Lambda }_i{\mathbf{V}}_i^x$$

$${\hat{A}}_f{\mathbf{V}}_q^x=$$

$$\frac{\mathrm{d}{\hat{A}}_f}{\mathrm{d}\sigma }{\mathbf{\kappa }}_r=-{\hat{A}}_f\frac{\mathrm{d}{\mathbf{\kappa }}_r}{\mathrm{d}\sigma }.$$

$${\hat{A}}_f{\mathbf{V}}_0^x=-\frac{\mathrm{d}{\hat{A}}_f}{\mathrm{d}\sigma }{\mathbf{U}}^x.$$

$${\hat{A}}_f{\mathbf{V}}_0^x={\hat{A}}_f\sum _r{a}_r\frac{\mathrm{d}{\mathbf{\kappa }}_r}{\mathrm{d}\sigma }$$

$${\mathbf{\kappa }}_0=\sum _r{a}_r\frac{\mathrm{d}{\mathbf{\kappa }}_r}{\mathrm{d}\sigma }$$

Full text

Modern models of excitability of cardiac cells include description of ionic channels in terms of deterministic Markov chains. The transition rates in these chains often vary by several orders of magnitude, which makes numerical simulations more difficult and necessitates development of specialized numerical approaches. We follow the usual wisdom that the small parameters in a mathematical model can be turned from an impediment into an advantage by development of asymptotic description exploiting those small parameters. The usual problem with experiment-derived (rather than postulated) models is that small parameters in them are not defined a priori but need to be identified, and sometimes this can be done equally plausibly in more than one way. In this paper, we show how the standard fast-slow asymptotic theory can be applied to the Markov chain models of ionic channels using one important model of this class as an example. We explore three selected ways of identifying small parameters in this model, and investigate the factors on which the utility of resulting asymptotics depends.

I Introduction

The bioelectricity is one of the driving forces of our life. Our mind and body is a manifestation of complex dynamics of electric impulses that carry information and trigger reactions in different organs of our body. The pulses of electrical excitation in heart are responsible for starting a chain of reactions, resulting in contraction of the cardiac muscle, causing the blood circulation. Understanding the detailed mechanisms of formation and propagation of electrical excitation can help in treatment and prevention of cardiac diseases.

The direct experimentation with living systems is difficult and rises many ethical issues, hence a mathematical description provides a valuable tool to gain insight and understanding of the internal working of the heart.

The elementary part of a cardiac excitation model are models of the ion-specific channel in the membrane, that close or open in response to the change in the transmembrane voltage. On the molecular scale, functioning of a single channel is an inherently stochastic process, which is adequately described as continuous time Markov chain (or Markov processes, as they are sometimes called). For most applications, however, the stochastic component is not essential, and it is sufficient to describe the behaviour of the channel in terms of the deterministic “master equation” for the probabilities of the channel to be in certain states, as functions of time. Simulation of resulting excitation models for single cells does not create problems; but when scaled to the tissue or whole-organ level, this becomes computationally expensive Plank et al. (2008). The computer technology is constantly improving, and recent publications describe results that would be unthinkable a few years before: e.g. Richards et al. (2013) simulated excitation in whole human heart with spatial resolution of 0.13 mm in nearly real time. However, that was done on a system with 1.6 million CPUs and a peak speed of 20 petaflops. For the moment, such systems are far from ubiquitous. It is therefore important to try and improve the computational methods for simulation of cardiac excitation models.

One significant factor of computational complexity is that the Markov chain models of ionic channels often involve processes on time scales differing by several orders of magnitude, i.e. are stiff. So a direct approach using explicit time steppers requires very small time steps, hence high computational demands. For instance, Bondarenko (2014), for a model involving a number of Markov chain channels, used time steps as short as 2 picoseconds—this is to be compared to the duration of the onset of an action potential of the order of 1 millisecond, and duration of the cardiac pulse of the order of 1 second. The natural alternative is generic implicit time steppers. However, approaches based on exploiting specific properties of cardiac excitation model present an attractive third possibility. In this paper, we explore one possible way to exploit specific properties of the Markov chain models of ionic channels. This is based on the traditional idea that a small parameter in the model can be turned from a hindrance into an advantage, by finding asymptotics in this parameter. For small parameters responsible for numerical stiffness, the adequate approach is singular perturbation theory, of fast-slow asymptotics.

The problem of stiffness of the description of the ion gates’ dynamics is not new for the Markov chain models, and was there already for the Hodkin-Huxley (HH) type “gate” models, starting from the Hodgkin and Huxley (1952). Out of many approaches to overcome stiffness in these models, arguably the most popular one follows the work by Rush and Larsen (1978). In some respects, it offers and “ideal” solution for the HH-type models, combining accuracy and stability. This approach exploits the fact that the equation controlling an HH-type gate is quasi-linear, in the sense that is is linear with respect to the gate variable, although it depends on other, “control”, variables in a non-linear way. Typically, the control variable is the transmembrane voltage, but in some channels it is concentration of Ca ions in addition or instead of the transmembrane voltage. If the control variables change negligibly during one time step, then “freezing” these variables allows one to write an “exact” solution. In the simplest formulation, the resulting computational scheme is first-order accurate, but the coefficient in the leading order term depends on the rate of the voltage (or whatever the control variable is) and does not depend on the stiffness of the equation with respect to the gate variable. This approach may be be formalized as fast-slow approach, where the control variable is slow and gate variable is fast. In that case, the leading-order solution for the gate variable, the “instant equilibrium”, corresponds to the limit when the time step is much longer than the characteristic time constant defined by the current values of the transition rate. However, the Rush-Larsen scheme in fact retains its accuracy when the time step is much less than, or is comparable to this characteristic time.

The Rush-Larsen approach has been so popular, it inevitably spawned a number of attempts to improve and/or extend it. For instance, Perego and Veneziani (2009) proposed how to increase its accuracy from first to second order, and Marsh, Ziaratgahi, and Spiteri (2012) suggested how to apply it to equations which are not quasi-linear, by linearising them for the duration of the time step. The Markov chain description poses a different sort of challenge: the equations are quasi-linear with respect to the Markov state occupancies so there is no need in linearization, but instead of one equation for a gate variable, there is a system of simultaneous linear equations. The “straightforward” Rush-Larsen-type approach, when each of the equations of the system is considered in turn, by freezing all other variables together with the control variables, has been proven effective for a number of examples (see Ref. [Plank et al., 2008] and references therein). A more radical approach for Markov chain generalization of the Rush-Larsen scheme, which utilized the exact solution for the whole linear system, thus avoiding extra errors caused by freezing some Markov states while updating others, was described in Refs. Stary and Biktashev, 2015a, b.

The methods discussed above offer efficient numerical schemes, but do not exploit the fact that the different processes within the same Markov chain may, and often do, have vastly different speeds. Taking these into account can, at least theoretically, offer further advantages. There have been a number of inspiring examples of this kind in literature. For instance, Hinch et al. (2004) and Plank et al. (2008) used the concept of “rapid equilibrium”, exploiting rapid transition rates between some of the states of a Markov chain, to effectively “merge” the closely connected states into one “combined” state, for a number of Markov chain models. This approach can be formalized as a leading-order asymptotic in the classical fast-slow perturbation theory descending from the works by Tikhonov (1952) and Fenichel (1979), for a particular form in which a small parameter ${\varepsilon}$ appears in the equations: as a factor $1/{\varepsilon}$ in front of the transition rates between two selected Markov states.

In the present paper, we seek to analyse this sort of asymptotics in more detail, following the full formalism, rather than immediately getting to the answer by following the rather obvious, but still only intuitive “rapid equilibrium” argument. The motivation for a more detailed analysis includes possibilities of generalization of the asymptotic approach to the fast/slow separation cases other than pairs of fast reciprocal transition rates, and getting higher-order terms in asymptotics. Moreover, an important theoretical question is whether stiff Markov chain formulations of ionic channels can always be well described by the standard singular perturbation theory. The intrigue here comes from the fact that the Hodgkin-Huxley description of some cardiac channels, specifically the fast sodium current channel, may be only partly described by the Tikhonov asymptotics, in regards of the the activation, ‘$m$’-gate; whereas the attempts to treat the inactivation, ‘$h$’-gate as fast or slow compared to the transmembrane voltage are ineffective, and at the very least fail to describe some essential qualitative features, such as the maximum of the action potential which is lower than the reversal potential for the sodium ions and which is different in a single cell in a propagating wave in tissue, to name the simplest example Biktashev (2002); Biktashev and Suckley (2004a); Biktasheva et al. (2006a); Simitev and Biktashev (2011). Hence the answer for the Markov chain model of the same channel is far from obvious a priori and requires investigation.

To address the theoretical question posed above, we have chosen the Markov chain model of the fast sodium current developed by Clancy and Rudy (2002). This is not the stiffest model of existing models of the kind, but it is one of the most popular ones. As our study was methodological rather than practical, it was important for our choice that the transition rates, and related characteristic times, in this model are varied in wide ranges which makes identification of small parameters a nontrivial issue. In other words, our aim was not to identify examples when the standard asymptotics can successfully treat the problem of stiffness (such examples do exist, see above), but rather analyse cases when the standard approach fails, as a necessary step towards developing more adequate, non-standard approaches. An extra motive for the choice of the fast Na current in this context was that the Hogdkin-Huxley description of this current is known to require non-Tikhonov asymptotics, as discussed above.

The structure of the paper is as follows. Section II introduces the notation and main principles of the singular perturbation approach we are using. Section III discusses amendments required of this approach with account of the specifics of Markov chain models. Section IV presents a formalization of the process of identification of small parameters in experiment-based models, which we call parametric embedding. Section V introduces the Markov chain model of the fast sodium current which we use to apply the singular perturbations. Section VI presents the main results, coming out of a few different parametric embeddings of this model. This is concluded by discussion in Section VII. We also present an Appendix containing technical material which is required for reproducing the main results but not for their understanding.

II General Theory for Dimensionality Reduction

The singular perturbation theory is well known in a variety of different formulations. We mostly follow the terminology and notation used e.g. in Refs. Biktashev, 2003; Suckley and Biktashev, 2003; Biktashev and Suckley, 2004b; Biktasheva et al., 2006b; Biktashev et al., 2008, adjusting where necessary for our present purposes.

We consider an autonomous system of ordinary differential equations

[TABLE]

where ${\mathbf{{u}}},{\mathbf{f}},{\mathbf{h}}\in{\mathbb{R}}^{n}$, and ${\varepsilon}$ is a small positive parameter. We assume existence of a stable ${m}$-dimensional manifold $\{{\mathbf{{U}}}\}$ of equilibria of the unperturbed system, ${\varepsilon}=0$, i.e. ${\mathbf{f}}({\mathbf{{U}}})={\mathbf{0}}$, where ${\mathbf{0}}$ stands for the null vector, with coordinates ${\mathbf{{a}}}\in{\mathbb{R}}^{m}$, $1\leq{m}<{n}$, and looking for solutions of the perturbed system, ${\varepsilon}>0$, in the form

[TABLE]

where the perturbation of the solution ${\mathbf{{b}}}\in{\mathbb{R}}^{n}$ is orthogonal to the manifold, in the sense that

[TABLE]

where the vectors ${\mathbf{{V}}_{{i}}}({\mathbf{{a}}})$ are right eigenvectors of a Jacobian matrix ${\hat{F}}({\mathbf{{U}}})={\partial}{\mathbf{f}}/{\partial}{\mathbf{{u}}}|_{{\mathbf{{u}}}={\mathbf{{U}}}}$,

[TABLE]

and the summation index ${\ell}$ runs through the stable eigenvalues, ${\mathop{\mathrm{Re}}\left({\Lambda_{{\ell}}}\right)}<0$, ${\ell}={m}+1,\dots,{n}$, skipping the zero eigenvalues, ${\Lambda_{{k}}}$, ${k}=1\,\dots,{m}$, corresponding to the directions tangent to the manifold. Table 1 summarises the meaning of these and other index conventions as used throughout the text, subject to a small amendment in the next section.

The right eigenvectors corresponding to zero eigenvalues ${\Lambda_{{k}}}=0$ are tangent to the invariant manifold and can be found as

[TABLE]

We substitute (2) and (3) into (1), expand the nonlinear functions into their Taylor series and separate the components using left eigenvectors ${\mathbf{W}_{{i}}^{T}}$ as projectors. The detailed derivation is presented in the Appendix A. The final result reads as the following system of ODEs

[TABLE]

where in the right-hand side of the first equation we have kept the leading order term, ${\mathbf{W}_{{k}}^{T}}{\mathbf{h}}({\mathbf{{U}}})$, and the first-order correction ${\mathcal{F}}_{i}$, which works out as

[TABLE]

where

[TABLE]

Equations (6a) and (6b) are coupled through higher-order terms, and to complete the reduction, we need to eliminate ${\mathbf{{b}}}$. For the solution of the manifold coordinates ${\mathbf{{a}}}$ up to ${\mathcal{O}({\varepsilon}^{N})}$ it is sufficient to find the correction term ${\mathbf{{b}}}$ up to ${\mathcal{O}({\varepsilon}^{N-1})}$. The leading order term for the correction ${\mathbf{{b}}}$ in terms of ${\mathbf{{a}}}$ can be found by solving (6b) using the integrating factor method. The solution also requires Taylor expansion of the integrating factor and of the non-homogeneous term. This leads to

[TABLE]

which is to be substituted into the first-order term into (6a), which then becomes a closed equation for ${\mathbf{{a}}}$.

III Dimensionality Reduction for Time-Inhomogeneous Markov Chains

The master equation for Markov chain models of ionic channels can be written in the form

[TABLE]

Entries in the vector of dynamical variables ${\mathbf{{x}}}\in{\mathbb{R}}^{n}$ represent the probabilities, that an ion channel resides in a particular state. Entries of the transition matrix ${\hat{A}}\in{\mathbb{R}}^{{n}\times{n}}$ describe the conditional probabilities of a channel in one given state to transit to another state per unit of time, i.e. transition rates. In reality, the matrix ${\hat{A}}$ depends on other dynamic variables of the model, e.g. the transmembrane voltage, which in turn are affected by the dynamics of the Markov chain; however this is not essential for the formalism we describe here and we assume that ${\hat{A}}$ is an explicit function of time, just for simplicity of notation. The sum of the entries in the vector of dynamical variables is equal to $1$, i.e. it is a stochastic vector. This implies that the sum of the entries in each column of the transition matrix ${\hat{A}}$ has to be equal to [math]. This is achieved as the entries on the diagonal of the transition matrix are a sum of the entries out of the diagonal for each column of the matrix. This property together with the fact that the non-diagonal elements are non-negative constitutes the definition of ${\hat{A}}$ as a left-stochastic matrix.

To use the theory described in Section II we have to take into account one simplifying fact and two complications. The simplifying fact is that the system (10) is linear. The complications are, firstly, that the theory described in the previous section applies to an autonomous system, but the Markov chain in (10) has an explicit time dependence of the transition matrix ${\hat{A}}({t})$. Secondly, the theory requires a small parameter, however the Markov chain models contain transition rates determined experimentally, and identifiction of any small parameters in such a case is a separate task, sometimes nontrivial.

The first complication is dealt with using autonomisation, which means that we introduce an additional dynamical variable ${\sigma}$ to represent time (henceforth referred to as “autonomous time”). Then the vector of dynamical variables is

[TABLE]

and the dynamic equation is

[TABLE]

Note that system (11) is no longer linear unless the function ${\hat{A}}({\sigma})$ is a constant.

To address the second complication, we introduce the small parameters artificially in an empirical procedure we call parametric embedding, which is discussed in detail in the next section. For now it is important that as a result, we can split the transition rates matrix ${\hat{A}}$ into a fast part ${\hat{A}_{f}}$ and a slow part ${\hat{A}_{s}}$, and the difference between them is identified by the small parameter ${\varepsilon}$ appearing as

[TABLE]

We restrict consideration to the embeddings in which ${\hat{A}_{f}}$ and ${\hat{A}_{s}}$ are left-stochastic matrices. We assume that the fast matrix ${\hat{A}_{f}}({\sigma})$ is diagonalizable, and introduce the eigenvalues ${\lambda_{{q}}}({\sigma})$ and the right eigenvectors ${\bm{\kappa}_{{q}}}({\sigma})$:

[TABLE]

(and drop from now on the dependence on ${\sigma}$, for brevity). We assume that for all ${\sigma}$, matrix ${\hat{M}_{1}}$ has a full set of eigenvectors, the first ${m}\geq 1$ of the eigenvalues are zero, and the remaining are all real (and of course negative)111Note that diagonalizability and reality of the eigenvalues of the full transition rate matrix ${\hat{A}}$ can be guaranteed under the assumption of detailed balance Anderson (1991), and ${\hat{A}_{f}}=\lim_{{\varepsilon}\to 0}\left({\varepsilon}{\hat{A}}\right)$. . Correspondingly, we introduce also the left eigenvectors ${\bm{\rho}_{{q}}}$,

[TABLE]

Differentiation of the last identity with respect to ${\sigma}$ yields a relationship that will be useful:

[TABLE]

We transform the system (11) into fast time ${\tau}={t}/{\varepsilon}$ to get a system

[TABLE]

This can be considered in the format of (1) with

[TABLE]

The dimensionality of the autonomized system (15) is ${n}+1$; we keep the upper value of the corresponding indices as ${n}$ but reserve the value ${0}$ for the time variable ${\sigma}$; this is where parameter ${\iota}$, designating the minimal value of the eigenvalues’ indices in  Table 1, becomes [math].

The manifold of equilibria in this case is in fact a linear subspace of ${\mathbb{R}}^{{n}+1}$ which is the hull of the one-dimensional subspace corresponding to the time coordinate ${\sigma}$ and the kernel of the fast matrix ${\hat{A}_{f}}$:

[TABLE]

where ${a}_{0}={\sigma}$, ${\mathbf{{a}}}={\begin{bmatrix}{a}_{1},\dots a_{m}\end{bmatrix}}{{}^{T}}$ and

[TABLE]

To construct the reduced system, we need to find the Jacobian of function ${\mathbf{f}}$ and solve the eigenvalue problem. The Jacobian is easily found as

[TABLE]

Let us denote the components of the eigenvectors as

[TABLE]

Substituting (19) and (20) into ${\Lambda_{{i}}}{\mathbf{{V}}_{{i}}}={\hat{F}}{\mathbf{{V}}_{{i}}}$, we get

[TABLE]

Let us consider separately the cases ${\mu}_{i}=0$ and ${\mu}_{i}\neq 0$. For ${\mu}_{i}=0$, equation (21b) becomes

[TABLE]

which is the definition of an eigenvalue problem for matrix ${\hat{A}_{f}}$, so we can take ${\mathbf{{V}}^{{x}}_{{q}}}={\bm{\kappa}_{{q}}}$ and ${\lambda_{{q}}}={\Lambda_{{q}}}$, for ${q}={1},\dots,{n}$, out of which the first ${m}$ are zero eigenvalues. For ${r}={1},\ldots,{m}$, we have ${\Lambda_{{r}}}=0$, and differentiation of (13) with respect to ${\sigma}$ gives

[TABLE]

We find one more eigenpair for the case ${\mu}_{0}\neq 0$. Let us normalise the corresponding eigenvector so that ${\mu}_{0}=1$. Then to satisfy equation (21a) we must have ${\Lambda_{{0}}}=0$, and (21b) becomes

[TABLE]

If we substitute (18) into (24) and use (23), we get

[TABLE]

so we can choose

[TABLE]

(this is choice is of course non-unique because the zero eigenvalue has multiplicity ${m}+1$).

The left eigenvectors are treated similarly. To summarise the results, the eigenvalues ${\Lambda_{{i}}}$ and eigenvectors ${\mathbf{{V}}_{{i}}}$, ${\mathbf{W}_{{i}}}$ of the Jacobian in the time-extended system are related to those ${\lambda_{{q}}}$, ${\bm{\kappa}_{{q}}}$, ${\bm{\rho}_{{q}}}$ of the transition rate matrix via the following relationships:

[TABLE]

With these the time-component of the “leading-order term” works out as ${\mathbf{W}_{{0}}^{T}}{\mathbf{h}}({\mathbf{{U}}})=1$, as should be expected. For the Markov chain subspace, we use (14), (16), (18) and (27) to get

[TABLE]

and then (9) gives the formula for the components of the transversal correction,

[TABLE]

and the transversal correction itself as

[TABLE]

Finally, the first-order accurate reduced system of ODEs is given by (6a), leading to

[TABLE]

This result can be written in the matrix form as

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and dash ′ stands for differentiation with respect to ${\sigma}$.

IV Parametric embedding

To address the second complication, we introduce the small parameters artificially in a procedure known as parametric embedding, previously introduced in Refs. Suckley and Biktashev, 2003; Biktasheva et al., 2006b; West et al., 2015. This procedure is a formalization of the replacement of a small constant with a small parameter.

Definition 1

We will call a system

[TABLE]

depending on parameter ${\varepsilon}$, a one-parametric embedding of a system

[TABLE]

if $f(u)\equiv F(u,1)$ for all $u\in{\mathop{\mathrm{dom}}\left(f\right)}$. If the limit ${\varepsilon}\to 0$ is concerned then we call it an asymptotic embedding.

The typical use of this procedure has the form of a replacement of a small constant with a small parameter. If a system contains a dimensionless constant $a$ which is “much smaller than 1”, then replacement of $a$ with ${\varepsilon}a$ constitutes a 1-parametric embedding; and then the limit ${\varepsilon}\rightarrow 0$ can be considered. In practice, constant $a$ would more often be replaced with parameter ${\varepsilon}$ rather than multiplied by it, but mathematically speaking, in the context of ${\varepsilon}\rightarrow 0$ and $a={\mathrm{const}}\neq 0$, these two ways are formally equivalent. This explains the paradoxical use of a zero limit for a parameter whose true value is one.

In some applications, the “small parameters” appear naturally and are readily identified. However, this is not always the case, and in complex systems identification of adequate small parameters may be a task in itself, which is where the formalization of this procedure can be helpful. In the context of the definition above, it is important to understand that there are infinitely many ways a given system can be parametrically embedded, as there are infinitely many ways to draw a curve $F(u;{\varepsilon})$ in the functional space given the only constraint that it passes through a given point, $F(u;1)=f(u)$. In terms of asymptotics, which of the embeddings is “better” depends on the qualitative features of the original systems that need to be represented, or classes of solutions that need to be approximated.

If a numerical solution of the system can be found easily, then there is a simple practical recipe: to look at the solutions of the embedding at different, progressively decreasing values of the artificial small parameter ${\varepsilon}$, and see when the features of interest will start to converge. If the convergent behaviour is satisfactorily similar to the original system with ${\varepsilon}=1$, the embedding is adequate for these features.

To summarize, we claim that identification of small parameters in a given mathematical model with experimentally measured functions and constants will, from the formal mathematical viewpoint, always be arbitrary, even though in the simplest cases the choice may be so natural that that this ambiguity is not even realized by the modeller, and that “validity” of such identification can be defined only empirically: if the asymptotics describe the required class of solutions sufficiently well. The exceptions may be when the asymptotic series are in fact convergent, the approximation errors can be estimated a priori, but this is rare.

In the subsequent text, slightly abusing the above definition for the sake of brevity, we refer as “embedding” to particular instances of one-parametric embedding of a given system for a selected value of the parameter ${\varepsilon}$. The overall structure of the embeddings is always like in (12), and the difference is in the choice of the matrices ${\hat{A}_{f}}$ and ${\hat{A}_{s}}$.

V Definition of the Markov Chain model of the fast sodium

current (INa)

We apply the asymptotic theory described above to a Markov chain model of the fast sodium current developed by Clancy and Rudy (2002) (we consider the wild-type version). The relevant part of the model has the form

[TABLE]

where ${V_{m}}$ is the transmembrane voltage, ${g_{\mathrm{Na}}}$ is the maximal conductance of the fast Na current, ${E_{\mathrm{Na}}}$ is the reversal potential of the Na ions due to the transmembrane difference in the concentration of these ions, ${O}$ is one of the components of the vector ${\mathbf{{x}}}$ representing the fraction of open channels, corresponding to the open state of the fast Na current channels, ${I}_{\ell}$ represent all other transmembrane currents, and the vector ${\mathbf{X}}$ comprises all other dynamic variables of the model, such as other ionic channels, concentrations etc.

Fig. 1(a) shows the diagram of the Markov chain. We find it more convenient to rename the dynamic variables, i.e. the names of the states of the Markov chain, as reported in Ref. [Stary and Biktashev, 2015a]: these are single-letter names, as opposed to the original names in Ref. [Clancy and Rudy, 2002] which use up to three symbols. The only state in the model that corresponds to the the channel being open is ${O}$, and this name coincides with the nomenclature used by Clancy and Rudy. So, for this Markov chain we have ${n}=9$ and the state vector

[TABLE]

According to the diagram of fig. 1(a), the transition rate matrix has the structure

[TABLE]

Here and elsewhere in transition rates matrices, in the interests of saving space, we do not show diagonal elements and replace them with ${*}$: they are uniquely defined by the condition that the sum of elements in each row should vanish. So e.g. the top left diagonal element in the above matrix is $-{\overset{\rightharpoonup}{PO}}-{\overset{\rightharpoonup}{UO}}$, and the bottom right element is $-{\overset{\rightharpoonup}{VW}}$.

All the transition rates in ${\hat{A}}$ are functions of the transmembrane voltage ${V_{m}}$. Their exact definitions can be found in the original publication Clancy and Rudy (2002) (see also Refs. Stary and Biktashev, 2015a; Starý, 2016) and we do not present them here; however fig. 1(b) gives a graphical illustration of the magnitudes of these rates in the physiological range of ${V_{m}}$. In that figure, we use the sum of the transition rates between two states as the measure of the speed of their connection, i.e. ${{\mathrm{i}}\cdot{\mathrm{j}}}{\;\triangleq\;}{\overset{\rightharpoonup}{{\mathrm{i}}{\mathrm{j}}}}+{\overset{\rightharpoonup}{{\mathrm{j}}{\mathrm{i}}}}$. Indeed, it is this quantity that determines the speed with which the dynamic equilibrium between the two states is reached if occupancies of all other states are fixed.

VI Embeddings of the INa Model

Fig. 1(b) allows one to see what transition rates may be considered “fast” and thus included into ${\hat{A}_{f}}$. For instance, connections $PO$, $QP$, $TU$ are relatively fast in the whole range of voltages, connection $VW$ is always slow, whereas connection $OU$ is fast at high ${V_{m}}$ but not so high at low ${V_{m}}$, and connection $VU$ is somewhat intermediate between the group of clearly fast connection and the group of clearly slow connections. In accordance with the above discussed formal definition and informal semantics of the concept of embedding, we intend to treat the question of which connections can or should be considered fast as strictly empirical, so Fig. 1(b) does not provide the ultimate answer to this question, but merely the possible directions of search. For simplicity, we always construct ${\hat{A}_{f}}$ by including into it reciprocal pairs of transition rates, and completing the diagonal elements to ensure ${\hat{A}_{f}}$ is left-stochastic. As a result, any nonzero non-diagonal element of ${\hat{A}}$ is always included into ${\hat{A}_{f}}$ and/or into ${\hat{A}_{s}}$, which guarantees that ${\hat{A}_{f}}$ and ${\hat{A}_{s}}$ are left-stochastic. We assess the quality of an embedding by how well it approximates the transients of the Markov states in a typical solution, a standard action potential; and of course of all the states the most important is the open state ${O}$.

We have tried a number of different combinations of reciprocal transition rates for ${\hat{A}_{f}}$. Not all such combinations pass the embedding test, i.e. give reasonable approximation of the original solution in the limit ${\varepsilon}\to 0$. In particular, the would-be “straightforward” solution to consider as “fast” all the transition rates that appear so in fig. 1(b), does not workStarý (2016). Figure 2 shows results of simulation of some of the more successful of those combinations. The ${I_{\mathrm{Na}}}$ model was extracted from the authors code Clancy and Rudy (2002). The simulation of the model were driven by recorded values of ${V_{m}}({t})$ during a standard action potential from a single-cell simulation. That means, we have performed a simulation of the original full model (35) once, and the resulting function ${V_{m}}({t})$ was then used for computations of the embedded version of only the subsystem  (35b), in which ${V_{m}}({t})$ was considered given and fixed; in other words, performed “virtual voltage clamp” experiments. The time step in the simulation of ${I_{\mathrm{Na}}}$ was ${\Delta_{t}}=1\;\mu\mathrm{s}$. The original model is shown with red lines, the transition rates embeddings are shown for a value of ${\varepsilon}=0.1$.

As can be deduced from the figure, although the transition rates included in the embeddings have roughly the same orders of magnitude, their expected effect on the accuracy of approximation of the ${O}$ transient by asymptotic methods is rather different: the ${O}{P}$ embedding is relatively poor, the ${R}{Q}$ is somewhat better, whereas embeddings involving transitions between ${S}$, ${T}$ and ${U}$, any pair or all three, promises very good accuracy: the corresponding graphs are indistinguishable in the plot resolution. Note that this assessment heavily depends on the special role of the ${O}$ state, and would be completely different if we were more interested in another Markov state. For instance, for the ${S}({t})$ transient, the OP embedding promises good accuracy, and STU embedding is very poor. Obviously, it matters how close are the embedded rates to the state in question.

VI.1 OP-embedding

In this section we develop an example of a particular embedding of the transition rates between the states ${O}$ and ${P}$, i.e. rates ${\overset{\rightharpoonup}{OP}}$ and ${\overset{\rightharpoonup}{PO}}$. As seen from the above discussion, the empirical evidence suggests that the asymptotics of this embedding is not likely to give a good approximation, so the purpose of this exercise is mainly didactic, to demonstrate in detail the application of the general theory, including the first-order correction, on a simple example.

In this embedding, the transition matrix ${\hat{A}}$ is split according to (12) into the matrix of the slow transition rates

[TABLE]

and the matrix of the fast transition rates

[TABLE]

For the dimensionality reduction we need to calculate the eigenvalues and eigenvectors of the fast matrix ${\hat{A}_{f}}$. This will result to a number of zero eigenvalues corresponding to the zero part of the matrix. There will be also at least one zero eigenvalue ${\lambda_{{1}}}=0$ corresponding to the Markov chain since

[TABLE]

In fact, we have ${\lambda_{{r}}}=0$ for ${r}=1,\dots,{m}$, where ${m}=8$, and just one non-zero eigenvalue ${\lambda_{9}}=-({\overset{\rightharpoonup}{PO}}+{\overset{\rightharpoonup}{OP}})$. The corresponding right eigenvectors are

[TABLE]

where ${\mathbf{e}_{{i}}}$ is the standard notation for the column-vector which has ${i}$-th component equal to one and all other components equal to zero, and the left eigenvectors are

[TABLE]

We note that the left eigenvector ${\mathbf{1}}$ asserted by the identity (39) is a linear combination of these, namely ${\mathbf{1}}=\sum_{{i}=1}^{8}{\bm{\rho}_{{i}}}$. The choice of normalization for ${\bm{\kappa}_{1}}$ and ${\bm{\rho}_{1}}$ is motivated by the ease of interpretation of the slow variable ${a}_{1}$, which will transpire shortly below.

Now we are ready to substitute the specifics of the selected embedding into the equation (30b) describing the reduced model. The left eigenvectors are constant for all ${r}$, so their derivatives are zero, and $\{{\ell}\}=\{9\}$. Then upon substituting (28) into (30b) we get

[TABLE]

The differential equation for ${a}_{r}$ for ${r}=3,4,5,7,8$ come out identical to the equations for the states ${R},{S},{T},{V},{W}$ from (36). This is because the first-order term vanishes as ${{\bm{\rho}_{{r}}}^{T}}{\hat{A}_{s}}({\sigma}){\bm{\kappa}_{9}}=0$ for these ${r}$. Hence we retain the same names for the corresponding components of the reduced model, as they had in the original model, and the vector of dynamic variables in the reduced system has the form

[TABLE]

where ${{\tilde{N}}}{\;\triangleq\;}{a}_{1}$, ${{\tilde{Q}}}{\;\triangleq\;}{a}_{2}$ and ${\tilde{U}}{\;\triangleq\;}{a}_{6}$.

The components ${r}=1,2,6$ in (42), that is differential equations for ${a}_{1}={{\tilde{N}}}$, ${a}_{2}={{\tilde{Q}}}$ and ${a}_{6}={\tilde{U}}$, will have nonzero first-order terms. According to ${\bm{\rho}_{1}}$ as given by (41), the new variable ${{\tilde{N}}}$ is just a sum of the old states occupancies ${O}$ and ${P}$; this is where the chosen normalization for ${\bm{\rho}_{1}}$ comes helpful. The names of the slow variables ${{\tilde{Q}}}$ and ${\tilde{U}}$ are motivated by the fact that according to (41) they map exactly to ${Q}$ and ${U}$ respectively, and the difference from the old variables is only in the first-order corrections in the reduced differential equations they obey.

Equation (18) then defines the relationship between the original and the reduced variables in the leading order, which in our case is

[TABLE]

where we define the fractions of the transition rates as

[TABLE]

We have only one stable eigenvalue in the present case, so equation (3) reduces to

[TABLE]

and equation (28), with account of ${\bm{\kappa}_{{r^{\prime}}}}^{\prime}=0$, ${r^{\prime}}\neq 1$, gives

[TABLE]

Then the leading-order transition matrix, according to (32,III) is

[TABLE]

where the new transition rates are defined as

[TABLE]

and the first-order correction to the transition matrix defined by (33,III) works out as

[TABLE]

The Markov chain of the ${I_{\mathrm{Na}}}$ channel is linked to the rest of the cell excitability model via the state ${O}$ which is the probability of the channel being open, so we need to compute ${O}$ in terms of the new dynamic variables. This is obtained from

[TABLE]

where ${{\mathbf{{U}}}^{x}}$ is given by (44) and ${{\mathbf{{v}}}^{x}}$ is given by (46), with (47) giving ${b}_{9}$. This leads to

[TABLE]

where

[TABLE]

Matrix ${\hat{M}_{1}}$ and coordinate ${b}_{9}$ depend on time derivatives of the transition rates, which in fact depend on the transmembrane voltage, hence the time derivative are to be calculated by the chain rule, e.g.

[TABLE]

Figure 3 shows the simulation results in the ${O}{P}$-embeddings and corresponding ${O}{P}$-reduction. The results should be compared against the original model shown by red lines. The simulations were done using extracted Markov chain model of ${I_{\mathrm{Na}}}$ driven by recordings of membrane voltage from whole cell simulations saved every $0.01\,\mathrm{ms}$ and interpolated as necessary. The state ${O}$ in reduced model was computed using (50,51). For comparison of the reduced model with the embedding, the occupancy of state ${{\tilde{N}}}$ in the embedded model was found as ${{\tilde{N}}}={O}+{P}$.

The simulations with leading order approximation (blue lines) show relatively large deviation from the original model. The first-order accurate asymptotic model computed for ${\varepsilon}=0.5$ (magenta lines) provides better approximation than only the leading order term, however the state ${O}$ in this approximation goes below zero, which does not make sense physically, as it represents a probability, so should be in the interval $[0,1]$: note that the generic asymptotic theory does not take into account these specifics.

The panel (d) shows the error norms computed using the following formula

[TABLE]

where ${\tilde{N}_{\mathrm{ref}}}$ is the reference solution obtained for a very small time step, and comparison is done for the interval of ${{t}_{\max}}=2$ ms of time-evolution. The error norms increase monotonically with ${\varepsilon}$ and show the convergence for the leading-order and first-order approximations as expected, which confirms the correctness of the formulas.

VI.2 STU-embedding and reduction of S, T and U into M

In this section we develop another approximation of the original system, which considers the transitions between states ${S}$, ${T}$ and ${U}$ as fast, which in asymptotics leads to their merger into a new state ${\tilde{M}}$. This choice is supported by the empirical embedding procedure as described in Section IV, details can be found in Ref. [Starý, 2016]. Now the matrix of fast transition rates is

[TABLE]

and the remaining, slow rates constitute the matrix

[TABLE]

The right eigenvectors corresponding to zero eigenvalue of this system can be chosen as

[TABLE]

The corresponding left eigenvectors are

[TABLE]

With account of these, we can keep the names of the original dynamic variables for all states except ${S}$, ${T}$, ${U}$, so the vector of states of the reduced system is

[TABLE]

These are all the ingredients needed for the derivation of the leading-order approximation. We have ${\bm{\rho}_{{i}}}^{\prime}=0$ for all ${i}=1,\dots,7$ so (32) gives the leading-order transition rate matrix for the reduced model as

[TABLE]

with the new transition rates defined as

[TABLE]

These expression use the notation ${\gamma_{\mathrm{ijkl}}}$ as an abbreviation for

[TABLE]

The original coordinates are recovered from the reduced one by

[TABLE]

As can be seen in fig. 4 below, the quality of the approximation obtained with these asymptotics, is very good. This was of course to be expected based on the results of the empirical embedding study, as discussed above.

VI.3 Embedding and reduction of R and Q states of STU-reduction into L

In this section, we investigate how one can build on the success of the ${S}{T}{U}$ embedding and achieve further reduction. As we have already considered the ${O}{P}$ reduction above, we now consider ${R}{Q}$ reduction. That is, we consider the transition rates between ${R}$ and ${Q}$ as fast, which will lead to the merger of these two states into a new state ${\tilde{L}}$. So in the context of the present section, the “original model” is defined by the matrix (58), which will now be called ${\hat{M}}$, while the fast matrix in the new embedding is

[TABLE]

and the slow matrix is

[TABLE]

Acting as before, we find the right eigenvectors of ${\hat{M}_{f}}$ corresponding to zero eigenvalue as

[TABLE]

the corresponding left eigenvectors as

[TABLE]

and we set the names of the components of the reduced vector as

[TABLE]

The resulting leading-order reduced matrix works out as

[TABLE]

with the new transition rates defined by

[TABLE]

where

[TABLE]

It is easily seen that the resulting reduced model (65) will be the same if, instead, we do the $QR$ reduction first and $STU$ reduction second, or do them simultaneously, i.e. include all of ${S\cdot T}$, ${T\cdot U}$ and ${Q\cdot R}$ in ${\hat{A}_{f}}$ in the first place.

Fig. 4 presents the results of the ${Q}{R}$-${S}{T}{U}$, together with the previously considered ${O}{P}$ and ${S}{T}{U}$ reductions. In these simulations, the ${I_{\mathrm{Na}}}$ channel model was not driven by the recorded ${V_{m}}({t})$ as before, but rather was part of the full cell model (35). The full original model and the three reduced version were run in the same protocol, which included stimulation with a period of one second, starting from ${t}=1\,\mathrm{ms}$ (this was done in order to be able to show the time in panel (a) in the logarithmic scale). One can see that the reduced models are indistinguishable from the full model except for the upstroke of the action potential. The upstroke of the fifth action potential is shown in detail in panels (b) and (c), for the probability of the ${I_{\mathrm{Na}}}$ channel being open, and the resulting value of this current. We see that the results generally agree with what could be expected from the empirical embedding studies illustrated in fig. 2. Namely, the ${O}{P}$ embedding gives a rather poor approximation, the ${Q}{R}$-${S}{T}{U}$ embedding is slightly better, while ${S}{T}{U}$ is very good.

VII Discussion

Fig. 5 summarises the Markov chain models occurring as a result of the three asymptotics we have considered: this is to be compared with the original scheme shown in fig. 1.

Asymptotic reduction based on time scale separation can pursue at least two different goals: reducing the number of dynamic equations, and reducing stiffness of those equations. The reductions considered in this paper are not particularly impressive in terms of reducing the number of equations: we have reduced by maximum of three out of nine, which is even less significant in comparison with the number of other dynamic equations in a typical model of an excitable cell, beyond the Markov chain of the ${I_{\mathrm{Na}}}$ channel. However, in practical applications the main goal is the other one: reducing the stiffness. To achieve a simple practical estimate of this characteristic, we measured the stiffness of the model by the maximum time step size ${\Delta_{t}}$ which provides a stable solution using the forward Euler solver for the isolated ${I_{\mathrm{Na}}}$ model driven by a recorded action potential. The original full model allows the time step of about ${\Delta_{t}}\approx 0.04\,\mathrm{ms}$ for stable computations; an increase above that leads to numerical instability. In comparison to that, all three models considered allow ${\Delta_{t}}\approx 0.044\,\mathrm{ms}$, i.e. a rather modest improvement. The limited progress in this is due to the fact that in all three examples considered, we have included in the embedding only some of the fastest transition rates. And even in these cases, we have seen that asymptotic removal of some of the fast processes affects the accuracy of computations. Even though these effects are seen only during the upstrokes of the action potential, these upstrokes are of principal significance as they determine the conduction velocity in spatially-distributed simulation, and therefore also the more delicate and more important phenomena such as conduction block, wavebreaks etc. Hence further increase of the number of the reduced degrees of freedom does not seem to be an answer. Further research is of course needed to establish that with certainty, but, as already noted above, e.g. the “straightforward” approach embedding all the transition rates that appear “fast” in fig. 1 does not yield a satisfactory approximation Starý (2016).

From a more theoretical viewpoint, stiffness can sometimes be characterized by the eigenvalues of the system; in particular, the upper limit of the integration step is mainly affected by the eigenvalue with the largest absolute value. In fig. 6 we plot those absolute values for the original model and the reduced models, as functions of the transmembrane voltage. We see that whereas $OP$ reduction somewhat reduces stiffness at the lower end of the ${V_{m}}$ scale, it has virtually no effect at the upper end. On the contrary, $QR$-$STU$ reduction noticeably reduces stiffness at the upper end, without changing it at the lower end.

So, as far as the question posed in the introduction is concerned, the results obtained here seem to suggest that Tikhonov asymptotic structure, which implies fixed distribution of the roles of “fast” and “slow” variables, or, in this particular class of applications, rather “fast” and “slow” transition rates, may not be quite adequate for the this particular model of the fast sodium current, and some non-Tikhonov parametric embedding may be more fruitful, say when transition rates are considered fast in one range of ${V_{m}}$ and slow in the complementary range, possibly with the asymmetry between the reciprocal rates taken into account.

An alternative approach, which has proved to be more practical than the one considered here, has been described in our previous works Stary and Biktashev (2015a, b), dubbed “exponential solvers”. However, that approach is purely numerical and does not explicitly take into account the fast-slow structure of the model, hence an asymptotic approach seems to have an a priori advantage, which ought to have been explored. We hope that the present study fills this gap to a certain extent.

An attractive possibility to improve the accuracy of the asymptotics and hence to open the way to further decrease the number of equations and reduce the stiffness, seems to be using higher-order asymptotics. We have explored this only in one of the three examples, but it already shows that (i) the algebraic complexity of the resulting formulas increases considerably, (ii) more significantly, some improvement in accuracy is devalued by the fact that the resulting model, unlike the leading-order asymptotics, no longer behaves as a “proper” Markov chain: the vector of dynamic variables is not guaranteed to remain stochastic, in particular, it can easily lead to negative values of the state occupancies. This happens because we have used the asymptotic theory which was designed for generic systems and is not tailored for the specific requirements of Markov chains. Hence another possible way for improvement may be in developing higher-order asymptotics strictly within the class of Markov chains.

Acknowledgements.

VNB gratefully acknowledges the current financial support of the EPSRC via grant EP/N014391/1 (UK) TS acknowledges financial support of the University of Exeter via PhD Studentship and of the EPSRC via grant EP/N024508/1.

Appendix A Derivation of Reduced System

We use the Taylor expansion for the functions ${\mathbf{f}}({\mathbf{{u}}})$ and ${\mathbf{h}}({\mathbf{{u}}})$, such that, after substitution of the sought solution (2), we get the expression on the right hand side of (1) as

[TABLE]

The first term ${\mathbf{f}}({\mathbf{{U}}})=0$ by assumption, and in the second term we note that the derivatives constitute the Jacobian matrix and expand the ${\mathbf{{v}}}$ according to (3). Then the previous expression (67) rewrites as

[TABLE]

We substitute the sought solution also to the left hand side of (1). We use the knowledge of eigenvectors corresponding to zero eigenvalues from (5) and expand the perturbed term (3). Then using a chain rule for the derivative of ${\mathbf{{U}}}({\mathbf{{a}}})$ and ${\mathbf{{V}}_{{\ell}}}({\mathbf{{a}}})$ we get

[TABLE]

Combining the right-hand side given by (67) and the left-hand side given by (69), we rewrite (1) as

[TABLE]

Multiplying the equation by the adjoint eigenvectors ${\mathbf{W}_{{i}}^{T}}$ gives

[TABLE]

Considering separately the zero and the stable eigenvalues then yields equations (6a) and (6b) respectively.