Abstract
We study Hopf bifurcation solutions to the Monodomain model equipped with FitzHugh-Nagumo cell dynamics. This reaction-diffusion system plays an important role in the field of electrocardiology as a tractable mathematical model of the electrical activity in the human heart. In our setting the (bounded) spatial domain consists of two subdomains: a collection of automatic cells surrounded by collections of normal cells. Thus, the cell model features a discontinuous coefficient. Analytical techniques are applied to approximate the time-periodic solution that arises at the Hopf bifurcation point. Accurate numerical experiments are employed to complement our findings.
1. Introduction
Several important properties of cardiac tissue can be adequately modeled by today’s detailed ionic models. These quantitative cell models consist generally of dozens of ODEs and do not easily lend themselves to mathematical analysis; for a comprehensive introduction to the field, see ,for example, the review paper [1]. There are, however, simple phenomenological models suitable for analysis of nonlinear phenomena, for example, the FitzHugh-Nagumo model [2] considered in this paper. These models capture the qualitative characteristics of electrical propagation in the heart. Moreover, they are both computationally and mathematically tractable and may be a competitive alternative to ionic models for many important problems in electrocardiology. In fact, phenomenological models of FitzHugh-Nagumo type have often successfully predicted the behavior of cardiac tissue [3].
The present paper concerns the FitzHugh-Nagumo model for pulse propagation in a continuum of heart cells, where the myocytes are diffusively coupled as described by the Monodomain model [4, 5]. We are interested in the interaction of a small collection of automatic (self-oscillatory) cells surrounded by finite populations of excitable but stable cells. The Sinoatrial (SA) node is a group of self-oscillatory cells that initiate each normal heart beat. It might happen that collections of cells other than the SA node (or the Atrioventricular node) become self-oscillatory. A collection of automatic cells that normally are not self-oscillatory is called an ectopic focus [6], and the appearance of it might result in life-threatening cardiac arrhythmias. To model this scenario a discontinuous coefficient will be incorporated into the cell kinetics in the same way as in the papers [7, 8]. These papers investigated the influence of oscillatory and excitable regions on wave propagation in the myocardial tissue; numerical computations for the FitzHugh-Nagumo model on a one-dimensional medium were performed to obtain results on the propagation characteristics. In another paper [9], both experiments and numerical simulations using an ionic model of cell kinetics were used to enlighten the role of automaticity, heterogeneity, and diffusive coupling for the generation of ectopic waves in myocardium.
Stability properties of equilibrium solutions to equations modeling an ectopic focus have been investigated in, for example, [6, 10]. In [6] the FitzHugh-Nagumo model was used to find the critical size of a pacemaker. The model was considered with spherical symmetry on an infinite domain; analysis revealed qualitative relationships between the stability of the steady state solution and the size of the automatic region, and the strength of the diffusive coupling. In a previous paper [11] we have determined quantitatively when the steady solution to the FitzHugh-Nagumo model in our setting loses stability; model parameters were the diffusion coefficient, the size of the automatic region and the strength of automaticity. More precisely, the critical diffusion coefficient in the system was identified; if the diffusion is smaller than this number, the equilibrium solution loses stability. A weaker coupling (less diffusion) between heart cells may be caused by fibrosis [12].
The purpose of the present paper is to pursue the research beyond the question of the influence of model parameters on the stability of the steady equilibrium solution, as studied in [6, 10, 11]. Here we want to take nonlinear effects into play as well and verify analytically the small oscillating solutions that are observed numerically when the steady state breaks. Furthermore, we set out to find analytical approximations of these new solutions to the model on a bounded domain. By doing so, we contribute with expressions that facilitate the study of the influence of model parameters on the behavior of the cardiac tissue in the presence of a self-oscillatory clump of cells. The transition of stable solutions occurs at critical points having purely imaginary eigenvalues, and thus we are concerned with Hopf bifurcation; see for example, [13]. The bifurcation parameter is a small perturbation of the critical diffusion coefficient. In [14, 15] Hopf bifurcations for simplified FitzHugh-Nagumo models for nerve conduction were analyzed. These papers were followed by the paper [16] considering the full FitzHugh model on an infinite domain. However, that paper concerned the model in a setting different from ours, but relevant for propagating pulses in nerve conduction: the FitzHugh parameters were independent of position, a stimulus current was applied, and the bifurcation parameter was a disturbance to the current strength. Analytical approximations of the solutions were not at hand for that problem setup.
This paper is outlined as follows. In the next section the preliminaries needed to calculate the bifurcating solution are worked out. Then, the steps to obtain the approximative solution are carried out. This is followed by a section with numerical experiments, and finally a short summary concludes the paper.
2. Preliminaries
We consider the Monodomain model together with the FitzHugh-Nagumo description of cell kinetics. Our purpose is to study a spatial domain where one region consists of automatic cells and another region has normal quiescent cells. The normal cells are surrounding the collection of automatic cells, which would act as pacemaker cells if they were isolated. The governing PDE-system models the transmembrane potential (the potential difference across the cell membranes) and a so-called recovery variable . The equations are here considered in the following form:
augmented with homogeneous Dirichlet boundary conditions on the closed domain , where relates to the size of the automatic domain and . Here and are dimensionless parameters arising from the derivation of the phenomenological FitzHugh-Nagumo cell model; see [2, 6] for the details. The diffusion coefficient comes from the Monodomain reduction of the Bidomain model, which is a more detailed description of heart tissue; compare [5]. The automaticity parameter depends on the position and describes automatic and normal cells:
where controls the automaticity of the automatic cells and is a parameter characterizing the normal cells.
In a previous paper [11] we have determined the diffusion coefficient for which stability is lost; this loss of stability might occur at a Hopf bifurcation point. More specifically, a number is found such that the eigenvalues of the linearized equilibrium system for (2.1) are purely imaginary, with . Next follows a short digression on the eigensystem of the linearized equilibrium equations; complre [11] for the details.
2.1. Eigensystem
Linearizing the equilibrium system around the null-solution gives the eigenvalue problem:
on with homogeneous Dirichlet conditions. For and it is found that the eigenvalues where Hopf bifurcation possibly occurs are with
Belonging to the eigenvalue we have the eigenfunction
and has the eigenfunction . Here is a real constant that will be determined by normalization later on, and is the function
where and . Finally, the diffusion coefficient yielding these purely imaginary eigenvalues is ; the point is determined by the equation
which can be solved numerically by, for example, Newton's method.
In general, we can find the eigenvalue for a given diffusion coefficient through the characteristic equation
where
2.2. Notational Definitions
We are interested in the solutions to (2.1) with , where is the bifurcation parameter. Our goal is to approximate the time-periodic solutions that are born when the null-solution loses stability. To this end we will make use of the general framework worked out in [13]. First, some notation needs to be defined.
We will write (2.1) as
Here and the derivatives are
where for example, denotes the first derivative of evaluated at and acting on . The problem (2.10) is set in the context of Hilbert spaces; two function spaces are introduced. The first one is the Hilbert space with the scalar product
for two-component functions . The second is the space of -periodic functions having the scalar product
with for all .
2.3. The Adjoint Problem
We will need the eigenfunctions for the adjoint operator to . Recall that
Straightforward calculations show that and
for homogeneous Dirichlet boundary conditions. If we consider the eigenvalue problem and use that we obtain the scalar eigenvalue problem
Similarly, the adjoint eigenvalue problem can be reduced to a scalar problem by using that . This results in the scalar eigenvalue problem above with conjugated . Thus, we have that the eigenvalues for the two eigenvalue problems are the conjugate of each other. It also holds that the first component of the eigenfunction for the adjoint eigenvalue problem is, up to a constant factor, identical to the one for . For it is found that the adjoint eigenfunction is
where is given in (2.6) and is an unspecified complex constant. The adjoint eigenfunction for is given by .
3. The Bifurcating Solution
We are ready to look for the time-periodic solutions to (2.1) that arise when the null-solution becomes unstable. To begin with, consider the linearized equation with . Note that if we take the eigenfunction in (2.5) and let , then and its conjugate solve and , respectively. Letting we are therefore led to define the operator in [13]:
such that , and any solution to can be written as a linear combination of and .
Our task in the following is to find a first-order accurate analytical approximation of the unknown bifurcating solution to (2.10) for some bifurcation parameter . Projections to the adjoint eigenspace will be utilized to single out the approximative solution.
3.1. Principal Eigenvalue and Bifurcation Parameter
The null-solution loses stability if the real part of the principal eigenvalue becomes positive when the bifurcation parameter changes, say grows, from zero. Note that , where is the derivative of the principal eigenvalue with respect to , evaluated at . In a given problem setup with a certain parameter set we should check that the real part of is positive, thereby indicating an eigenvalue with positive real part for . Two different ways to calculate will be provided. Since the corresponding values obtained for a given bifurcation problem must be identical, equating the two formulas gives a condition that should be verified in order to increase the confidence in the calculations.
First, to obtain we can use the characteristic equation (2.8) for and calculate
and then evaluate this expression at to find . By using Maple an expression for is found; denoting it by it reads
Second, an alternative way to calculate is given in [13]. Consider the spectral problem , where at the function is identical to the eigenfunction given in (2.5). The equation resulting from taking the derivative of the spectral problem with respect to at is solvable if and only if
We turn to find the right-hand side in (3.4). The eigenfunctions and occurring in (3.4) need to be normalized such that . Let us start out by normalizing the eigenfunction given in (2.5) We could simply put in (2.5), but choose to determine it such that :
and so the sought constant can be written as
where we have defined
Next, the complex constant for the adjoint eigenfunction in (2.17) is determined from the condition . Using that we calculate
and find
Now, and
where is the integral
So, the condition arising from (3.3) and (3.4) may be written as
3.2. Series Expansion
We attempt to find the solution to (2.10) as a power series in a small parameter :
Here is the approximative solution, , is the bifurcation parameter, and is the time frequency.
A Fredholm [17] solubility condition is given in [13], First, from we may choose . That is,
It also holds that and
where . However, in order to solve this equation for and we have first to solve
for , which is not easily done. Since we are looking for a first-order approximation, we can make use of an alternative equation instead. This equation arises in the process of solving a complex-valued amplitude equation for the problem (2.10); it ignores third-order terms and reads [13]
where is given by the right-hand side of (3.12) and
We modify the equation to conform to the current series expansion and include the missing cubic term as well:
and here the cubic term reads
Recalling that and the coefficients in (3.19) may now be calculated. It is found that and
where we have defined the integrals
Now, solving (3.19) for and we obtain the bifurcation coefficient
The frequency coefficient is not significant for first-order accuracy but is given here for completeness:
3.3. Summary and Discussion of Analytical Results
The bifurcation parameter for a small parameter is given by
and the corresponding diffusion coefficient is with determined from (2.7). Here the frequency coefficient is and the integrals are given by
with in (2.2) and the function given in (2.6):
with and , In passing we note that the integrals can be conveniently calculated; the integrands are even and combinations of elementary functions. To first order in the bifurcating solution to (2.1) reads
where .
The eigenvalues of the linearized equilibrium operator evaluated at the null-solution can be obtained by solving the characteristic equation (2.8) and we found that they may be approximated by
and . A positive real part of and a positive bifurcation parameter mean that the steady state is unstable.
We remark that the model parameters are coupled and that a parameter study at the transition of stable solutions involves solving (2.7) to obtain the critical diffusion coefficient . Doing so, we observe, for example, that the spatial extension of the wave (3.28) into the surrounding tissue decreases with increasing . Thus, the normal cells become less excitable as increases. Likewise, if is decreased, the automatic cells become less self-oscillatory and a smaller portion of the normal tissue will excite. These observations can be explained from model parameters in the following way.
In the formal calculations below it is assumed that real part of is positive, that the imaginary part of the eigenvalue is nonzero, and that the small bifurcation parameter is positive. Then the steady state will bifurcate into the time-periodic solution approximated by (3.28). Without loss of generality we let and consider (2.7) with . Now, put and such that (2.7) can be written as
First, we consider the case where the automatic cells are not very self-oscillatory, that is, and is small. We assume a sufficiently large domain; that is, is big enough, such that we may approximate (3.30) with
Expanding the -function yields the equation
on ignoring higher-order terms. Solving for the diffusion coefficient we get , which is small since is small. Consider now the function above describing the spatial shape of the bifurcating solution. Note that is close to and that is large since is small. Since is close to the function will be small at the boundary between the collections of different cells, and because is large it will approach zero rapidly as the distance to the automatic collection grows. Thus, the wave will not propagate far into the surrounding tissue.
Next, consider the case where is large so that is big as well. We again assume that is sufficiently big; analogous to the previous case we then obtain equation (3.32) and . For large, and as before we have that is close to and that is large. The conclusion is that the surrounding tissue is not very excitable for large ; the greatly stable cells prevent the spreading of oscillations.
By using (2.7) we find that the critical diffusion coefficient grows when is increased or is decreased. In these cases will be small, and the bifurcating solution will approach zero more slowly as the distance to the automatic collection increases; therefore a larger part of the surrounding tissue will excite. So, the surrounding tissue is more excitable for smaller and the self-oscillatory cells have a stronger automaticity for larger ; the transition of solutions occurs at a stronger diffusive coupling, meaning that the steady solution requires a larger amount of diffusion in order to be stable for these cases.
A Matlab-script has been written to illustrate the propagation of the oscillations generated by the clump of automatic cells, into the surrounding collection of normal cells. Different values of and are considered. The output graphed in Figure 1 supports the above analysis. Observe that the propagation of oscillations into the surrounding group of normal cells can be substantial. For the largest value of and the smallest value of considered, the strength of the oscillations at the interface between the regions of different cells is around % of the amplitude of the wave. We remark that the spatial shape of the approximation given in (3.28) of the bifurcating solution to (2.1) is independent of the size of the bifurcation parameter in (3.25); the spatial shape is approximated by in (2.6) which is independent of the bifurcation parameter . However, the amplitude of the solution decreases as the bifurcation parameter approaches zero. As we let go to zero, goes to zero and approaches the null-solution—which is reasonable.
(a)
(b)
(c)
(d)
(e)
(f)
4. Numerical Experiments
The purpose of this section is to compare the analytical approximation (3.28) of the bifurcating solution to (2.1) to accurate numerical approximations of the solution.
Equation (2.1) with the parameters , , , and diffusion coefficient is considered on with and . From (2.7) we obtain that for purely imaginary eigenvalues , .
First, we check that condition (3.12) holds. The characteristic equation (2.8) enabled us to find an expression for . This expression is given in (3.3) and evaluating it for the current parameters gives . Denoting by the right-hand side of (3.12) and checking the condition shows that . Now, should agree to first order with obtained from (2.8) with . Note that since the real part of is positive, the null-solution is unstable provided that . For different we calculate given in (3.25) and check the errors . By using Müller's method, see [18], is found to an accuracy of around . In Table 2 the convergence results are listed. We remark that is positive and that the errors are of second order in as expected.
Next, we will compute a numerical solution to (2.1) and do a comparison to the analytical approximation of the solution. For the numerical approximation we discretize the interval , even ( in the present case), with points, such that and are centered between grid points. In this way the influence of the singularities at is diminished. It is readily found that we may choose , to accomplish that are located in the center of grid points. By discretizing (2.1) in space with homogeneous Dirichlet conditions we can formulate a semidiscrete system . Here is the numerical solution vector and is the -matrix of the classical second-order finite difference approximation of the second derivative acting on odd entries of (that is on ). The function reads
where is given in (2.1) and . We take points and integrate the semidiscrete system in time up to by the built-in Matlab solver ode15s with tolerance . The large time interval is chosen to assure that the numerical solution is converged; the amplitude and time-period of the solution is checked. A small initial condition is used in the numerical computation:
We take and obtain from Table 1 that and that the corresponding eigenvalue is . Therefore, the null-solution is unstable and will bifurcate into a stable time-periodic solution. The analytical approximation of this solution is (3.28) with . Figure 2 depicts the approximative solution (3.28) and the numerical approximation in six snapshots covering of the time-period . We observe a good agreement between the approximative solutions.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
In the final experiment we test the convergence of the analytical approximation towards the highly resolved numerical reference solution. The first four cases in Table 1 are considered. All numerical computations are run until the solution is deemed to have converged. The errors are measured over one full time period in the norm induced by the scalar product
where we recall that . The approximation (3.28) is translated in time to be in phase with the numerical solution. The results are recorded in Table 2; we conclude that the errors are of second order in and that the approximation (3.28) has the expected accuracy.
5. Summary
We have studied bifurcation solutions to the Monodomain model in connection with FitzHugh-Nagumo cell kinetics. By incorporating a discontinuous coefficient in this reaction-diffusion system we modeled the interaction of normal and automatic heart cells—the relevant scenario for the study of a so-called ectopic focus that might appear in the cardiac muscle.
The electrical activity of the heart is essential for its function. Waves generated by an ectopic focus might disturb this synchronized electrical flow and lead to lethal cardiac arrhythmias. In our model such an ectopic wave is realized as a disturbance of the resting state. The results show that the presence of a group of strongly automatic cells indeed stimulates the surrounding normal cells. The resting state is broken by a wave generated by the automatic cell collection and oscillations propagate into the normal heart tissue. Stronger automatic cells will trigger a larger portion of the surrounding myocytes.
The reduction of the reaction-diffusion system to algebraic equations valid around the critical point made the model more transparent and enabled efficient parameter studies to assess the behavior of the cardiac tissue. The influence of model parameters and the interplay between them were discussed. With sufficiently small diffusive coupling between myocytes the collection of automatic cells caused a breaking of the steady equilibrium solution into a small time-periodic wave. The shape, amplitude, and time-dependence of this solution as a function of a bifurcation parameter were found in closed form to an adequate accuracy. The output from numerical computations agreed very well with our theoretical results.