# Bifurcating Solutions to the Monodomain Model Equipped with FitzHugh-Nagumo Kinetics

Authors: Robert Artebrant
Journal: Journal of Applied Mathematics
 (2009)

---

## 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.


---

## Body

## 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 potentialv (the potential difference across the cell membranes) and a so-called recovery variable w. The equations are here considered in the following form:(2.1)vt=δvxx+1ϵf(v,w,p(x)),wt=v-γw,f(v,w,p(x))=v(v+p(x))(1-v)-w,
augmented with homogeneous Dirichlet boundary conditions on the closed domain [-Ma,Ma], where a relates to the size of the automatic domain and M>1. Here ϵ>0 and γ>0 are dimensionless parameters arising from the derivation of the phenomenological FitzHugh-Nagumo cell model; see [2, 6] for the details. The diffusion coefficient δ>0 comes from the Monodomain reduction of the Bidomain model, which is a more detailed description of heart tissue; compare [5]. The automaticity parameter p depends on the position x and describes automatic and normal cells:(2.2)p(x)={α,|x|≤a,-β,a<|x|≤Ma,
where α>0 controls the automaticity of the automatic cells and β>0 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 η0=δ/a2 is found such that the eigenvalues of the linearized equilibrium system for (2.1) are purely imaginary, λ±=±iω0 with ω0∈ℝ. 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:(2.3)δv′′+1ϵ(p(x)v-w)=λv,v-γw=λw
on [-Ma,Ma] with homogeneous Dirichlet conditions. For 0<ϵ<1/γ2 and α>ϵγ it is found that the eigenvalues where Hopf bifurcation possibly occurs are λ±=±iω0 with(2.4)ω0=1/ϵ-γ2.Belonging to the eigenvalueλ+ we have the eigenfunction(2.5)ϕ(x)=c·ξ(x)(1,1γ+λ+)=c·ξ(x)(1,ϵ(γ-iω0)),
and λ- has the eigenfunction ϕ¯(x). Here c is a real constant that will be determined by normalization later on, and ξ is the function(2.6)ξ(x)={coskx,|x|≤a,coskasinh(κa(M-1))sinh(κ(Ma-|x|)),a<|x|≤Ma,
where k=(α/ϵ-γ)/δ0 and κ=(β/ϵ+γ)/δ0. Finally, the diffusion coefficient yielding these purely imaginary eigenvalues is δ0=η0a2; the point η0 is determined by the equation(2.7)α/ϵ-γη0-arctan(β/ϵ+γα/ϵ-γ1tanh((M-1)β/ϵ+γη0))=0,
which can be solved numerically by, for example, Newton's method.In general, we can find the eigenvalueλ for a given diffusion coefficient δ=ηa2 through the characteristic equation(2.8)G:=gA(λ)η-arctan(gN(λ)gA(λ)1tanh((M-1)gN(λ)η))=0,
where(2.9)gA(λ)=-λ-1ϵ(λ+γ)+αϵ,gN(λ)=λ+1ϵ(λ+γ)+βϵ.


### 2.2. Notational Definitions

We are interested in the solutions to (2.1) with δ=(η0-μ)a2, where μ>0 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(2.10)ut=F(μ,u)=Fu(0∣u)+μFuμ(0∣u)+12Fuu(0∣u∣u)+16Fuuu(0∣u∣u∣u).
Here u=(v,w) and the derivatives are(2.11)Fu(0∣u)=(a2η0vxx+p(x)v-wϵ,v-γw),Fuμ(0∣u)=-(a2vxx,0),Fuu(0∣u∣u)=2(1-p(x)ϵv2,0),Fuuu(0∣u∣u∣u)=-6(v3ϵ,0),
where for example, Fu(μ∣u) denotes the first derivative of F evaluated at μ and acting on u. The problem (2.10) is set in the context of Hilbert spaces; two function spaces are introduced. The first one is the Hilbert space H=(L2(-Ma,Ma))2 with the scalar product(2.12)〈U,V〉=∫-MaMaU1V¯1+U2V¯2dx
for two-component functions U,V∈H. The second is the space of 2π-periodic functions having the scalar product(2.13)[U,V]=12π∫02π〈U,V〉(s)ds
with U(s),V(s)∈H for all s∈(0,2π).


### 2.3. The Adjoint Problem

We will need the eigenfunctions for the adjoint operatorL* to L=Fu(0∣·). Recall that(2.14)Lu=(a2η0vxx+p(x)v-wϵ,v-γw).
Straightforward calculations show that 〈Lu,z〉=〈u,L*z〉 and(2.15)L*u=(a2η0vxx+p(x)v+ϵwϵ,-v+ϵγwϵ)
for homogeneous Dirichlet boundary conditions. If we consider the eigenvalue problem Lu=λu and use that w=v/(λ+γ), we obtain the scalar eigenvalue problem(2.16)a2η0vxx+(p(x)ϵ-1ϵ(λ+γ)-λ)v=0.
Similarly, the adjoint eigenvalue problem L*u=λ¯u can be reduced to a scalar problem by using that w=-1/(ϵ(λ¯+γ)). 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 L. For λ¯+=λ-=-iω0 it is found that the adjoint eigenfunction is(2.17)ϕ*(x)=d·ξ(x)(1,-1ϵ(λ¯++γ))=d·ξ(x)(1,-γ-iω0),
where ξ is given in (2.6) and d is an unspecified complex constant. The adjoint eigenfunction for λ¯-=iω0 is given by ϕ¯*.


## 2.1. Eigensystem

Linearizing the equilibrium system around the null-solution gives the eigenvalue problem:(2.3)δv′′+1ϵ(p(x)v-w)=λv,v-γw=λw
on [-Ma,Ma] with homogeneous Dirichlet conditions. For 0<ϵ<1/γ2 and α>ϵγ it is found that the eigenvalues where Hopf bifurcation possibly occurs are λ±=±iω0 with(2.4)ω0=1/ϵ-γ2.Belonging to the eigenvalueλ+ we have the eigenfunction(2.5)ϕ(x)=c·ξ(x)(1,1γ+λ+)=c·ξ(x)(1,ϵ(γ-iω0)),
and λ- has the eigenfunction ϕ¯(x). Here c is a real constant that will be determined by normalization later on, and ξ is the function(2.6)ξ(x)={coskx,|x|≤a,coskasinh(κa(M-1))sinh(κ(Ma-|x|)),a<|x|≤Ma,
where k=(α/ϵ-γ)/δ0 and κ=(β/ϵ+γ)/δ0. Finally, the diffusion coefficient yielding these purely imaginary eigenvalues is δ0=η0a2; the point η0 is determined by the equation(2.7)α/ϵ-γη0-arctan(β/ϵ+γα/ϵ-γ1tanh((M-1)β/ϵ+γη0))=0,
which can be solved numerically by, for example, Newton's method.In general, we can find the eigenvalueλ for a given diffusion coefficient δ=ηa2 through the characteristic equation(2.8)G:=gA(λ)η-arctan(gN(λ)gA(λ)1tanh((M-1)gN(λ)η))=0,
where(2.9)gA(λ)=-λ-1ϵ(λ+γ)+αϵ,gN(λ)=λ+1ϵ(λ+γ)+βϵ.


## 2.2. Notational Definitions

We are interested in the solutions to (2.1) with δ=(η0-μ)a2, where μ>0 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(2.10)ut=F(μ,u)=Fu(0∣u)+μFuμ(0∣u)+12Fuu(0∣u∣u)+16Fuuu(0∣u∣u∣u).
Here u=(v,w) and the derivatives are(2.11)Fu(0∣u)=(a2η0vxx+p(x)v-wϵ,v-γw),Fuμ(0∣u)=-(a2vxx,0),Fuu(0∣u∣u)=2(1-p(x)ϵv2,0),Fuuu(0∣u∣u∣u)=-6(v3ϵ,0),
where for example, Fu(μ∣u) denotes the first derivative of F evaluated at μ and acting on u. The problem (2.10) is set in the context of Hilbert spaces; two function spaces are introduced. The first one is the Hilbert space H=(L2(-Ma,Ma))2 with the scalar product(2.12)〈U,V〉=∫-MaMaU1V¯1+U2V¯2dx
for two-component functions U,V∈H. The second is the space of 2π-periodic functions having the scalar product(2.13)[U,V]=12π∫02π〈U,V〉(s)ds
with U(s),V(s)∈H for all s∈(0,2π).


## 2.3. The Adjoint Problem

We will need the eigenfunctions for the adjoint operatorL* to L=Fu(0∣·). Recall that(2.14)Lu=(a2η0vxx+p(x)v-wϵ,v-γw).
Straightforward calculations show that 〈Lu,z〉=〈u,L*z〉 and(2.15)L*u=(a2η0vxx+p(x)v+ϵwϵ,-v+ϵγwϵ)
for homogeneous Dirichlet boundary conditions. If we consider the eigenvalue problem Lu=λu and use that w=v/(λ+γ), we obtain the scalar eigenvalue problem(2.16)a2η0vxx+(p(x)ϵ-1ϵ(λ+γ)-λ)v=0.
Similarly, the adjoint eigenvalue problem L*u=λ¯u can be reduced to a scalar problem by using that w=-1/(ϵ(λ¯+γ)). 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 L. For λ¯+=λ-=-iω0 it is found that the adjoint eigenfunction is(2.17)ϕ*(x)=d·ξ(x)(1,-1ϵ(λ¯++γ))=d·ξ(x)(1,-γ-iω0),
where ξ is given in (2.6) and d is an unspecified complex constant. The adjoint eigenfunction for λ¯-=iω0 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 ut=Lu with L=Fu(0∣·). Note that if we take the eigenfunction ϕ in (2.5) and let z=eiω0tϕ, then z and its conjugate solve zt=Lz and z¯t=Lz¯, respectively. Letting s=ω0t we are therefore led to define the operator in [13]:(3.1)𝕁0=-ω0dds+Fu(0∣·),
such that 𝕁0z=𝕁0z¯=0, and any solution u to 𝕁0u=0 can be written as a linear combination of z and z¯.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 μ>0. 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 λ+≈iω0+λμμ, where λμ is the derivative of the principal eigenvalue with respect to μ, evaluated at μ=0. 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 μ>0. 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 η=η0-μ and calculate(3.2)0=ddμG=Gλλμ+Gμ
and then evaluate this expression at μ=0 to find λμ=-Gμ/Gλ. By using Maple an expression for λμ is found; denoting it by λμG it reads(3.3)λμG=κa3k2[(κ2+k2)sinh2(κa(M-1))+Mκ2](κ2+k2)sinh(2κa(M-1))+2aκ[(k2+κ2)cosh2(κa(M-1))-k2M](1-iγω0).Second, an alternative way to calculateλμ is given in [13]. Consider the spectral problem λζ=Fu(μ∣ζ), where at μ=0 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 μ=0 is solvable if and only if(3.4)λμ=〈Fuμ(0∣ϕ),ϕ*〉,〈ϕ,ϕ*〉=1.We turn to find the right-hand side〈Fuμ(0∣ϕ),ϕ*〉 in (3.4). The eigenfunctions ϕ and ϕ* occurring in (3.4) need to be normalized such that 〈ϕ,ϕ*〉=1. Let us start out by normalizing the eigenfunction ϕ given in (2.5) We could simply put c=1 in (2.5), but choose to determine it such that 〈ϕ,ϕ〉=1:(3.5)〈ϕ,ϕ〉=∫-MaMaϕ1ϕ¯1+ϕ2ϕ¯2dx=c2(1+ϵ)∫-MaMaξ2(x)dx,
and so the sought constant can be written as(3.6)c=1(1+ϵ)I0,
where we have defined(3.7)I0:=∫-MaMaξ2(x)dx.
Next, the complex constant d for the adjoint eigenfunction ϕ* in (2.17) is determined from the condition 〈ϕ,ϕ*〉=1. Using that ω0=1/ϵ-γ2 we calculate(3.8)〈ϕ,ϕ*〉=∫-MaMaϕ1ϕ¯1*+ϕ2ϕ¯2*dx=2cd¯ϵω0(ω0+iγ)I0
and find(3.9)d=12cI0(1+iγω0).
Now, Fuμ(0∣ϕ)=-(a2ϕ1,xx,0)=-(a2cξ′′,0) and(3.10)〈Fuμ(0∣ϕ),ϕ*〉=-∫-MaMaa2cξ′′(x)d¯ξ(x)dx=-cd¯a2I1=-a2I12I0(1-iγω0),
where I1 is the integral(3.11)I1:=∫-MaMaξ(x)ξ′′(x)dx.
So, the condition arising from (3.3) and (3.4) may be written as(3.12)λμG=-a2I12I0(1-iγω0).


### 3.2. Series Expansion

We attempt to find the solution to (2.10) as a power series in a small parameter ϵ̂:(3.13)(uϵ̂(x,s,ϵ̂)μ(ϵ̂)ω(ϵ̂)-ω0)=∑n=1∞ϵ̂nn!(un(x,s)μnωn).
Here uϵ̂ is the approximative solution, s=ω0t, μ is the bifurcation parameter, and ω is the time frequency.A Fredholm [17] solubility condition is given in [13], First, from 𝕁0u1=0 we may choose u1=z+z¯=eisϕ+e-isϕ¯. That is,(3.14)u1(x,s)=2cξ(x)(cos(s),ϵ(γcos(s)+ω0sin(s))).
It also holds that μ1=ω1=0 and(3.15)3(-iω2+μ2λμ)+3[Fuu(0∣u1∣u2),z*]+[Fuuu(0∣u1∣u1∣u1),z*]=0,
where z*=eisϕ*. However, in order to solve this equation for μ2 and ω2 we have first to solve(3.16)𝕁0u2=-Fuu(0∣u1∣u1),[u2,z*]=0
for u2, 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](3.17)iω2-μ2λμ=-2σ0σ1+4σ1σ¯1+2σ2σ¯2/3iω0,
where λμ is given by the right-hand side of (3.12) and(3.18)σ0=12〈Fuu(0∣ϕ∣ϕ),ϕ*〉,σ1=12〈Fuu(0∣ϕ∣ϕ¯),ϕ*〉,σ2=12〈Fuu(0∣ϕ¯∣ϕ¯),ϕ*〉.
We modify the equation to conform to the current series expansion and include the missing cubic term as well:(3.19)iω2-μ2λμ2=-2σ0σ1+4σ1σ¯1+2σ2σ¯2/3iω0+3ρ,
and here the cubic term reads(3.20)ρ=16〈Fuuu(0∣ϕ∣ϕ∣ϕ¯),ϕ*〉.
Recalling that Fuu(0∣ϕ∣ϕ)=2/ϵ((1-p(x))ϕ12,0) and Fuuu(0ϕϕϕ)=-6/ϵ(ϕ13,0), the coefficients in (3.19) may now be calculated. It is found that σ2=σ1=σ0 and(3.21)σ0=c2d¯ϵI2=I22I0ϵ(1+ϵ)I0(1-iγω0),ρ=-c3d¯ϵI3=-I32I02ϵ(1+ϵ)(1-iγω0),
where we have defined the integrals(3.22)I2:=∫-MaMa(1-p(x))ξ3(x)dx,I3:=∫-MaMaξ4(x)dx.
Now, solving (3.19) for μ2 and ω2 we obtain the bifurcation coefficient(3.23)μ2=22γI22-3ϵω02I0I3(ϵaω0I0)2(1+ϵ)I1.
The frequency coefficient ω2 is not significant for first-order accuracy but is given here for completeness:(3.24)ω2=-43I22(ω02+γ2)ω03(1+ϵ)I03ϵ2=-43I22(ω0I0ϵ)3(1+ϵ).


### 3.3. Summary and Discussion of Analytical Results

The bifurcation parameter for a small parameterϵ̂ is given by(3.25)μ=ϵ̂22γI22-3ϵω02I0I3(ϵaω0I0)2(1+ϵ)I1,
and the corresponding diffusion coefficient is δ=(η0-μ)a2 with η0 determined from (2.7). Here the frequency coefficient is ω0=1/ϵ-γ2 and the integrals are given by(3.26)I0=∫-MaMaξ2(x)dx,I1=∫-MaMaξ(x)ξ′′(x)dx,I2=∫-MaMa(1-p(x))ξ3(x)dx,I3=∫-MaMaξ4(x)dx
with p in (2.2) and the function ξ given in (2.6):(3.27)ξ(x)={coskx,|x|≤a,coskasinh(κa(M-1))sinh(κ(Ma-|x|)),a<|x|≤Ma
with k=(α/ϵ-γ)/δ0 and κ=(β/ϵ+γ)/δ0, δ0=η0a2. 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(3.28)uϵ̂(x,t)=2ϵ̂cξ(x)(cos(ω0t),ϵ(γcos(ω0t)+ω0sin(ω0t))),
where c=1/(1+ϵ)I0.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(3.29)λ+≈iω0+λμμ,λμ=-a2I12I0(1-iγω0),
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 a=1 and consider (2.7) with δ=η0. Now, put σ=α/ϵ-γ and τ=β/ϵ+γ such that (2.7) can be written as(3.30)σδ-arctan(τσ1tanh((M-1)τδ))=0.
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, M is big enough, such that we may approximate (3.30) with(3.31)σδ-arctan(τσ)=0.
Expanding the arctan-function yields the equation(3.32)σδ-(π2-στ)=0
on ignoring higher-order terms. Solving for the diffusion coefficient we get δ=σ/(π/2-σ/τ)2, which is small since σ is small. Consider now the function ξ above describing the spatial shape of the bifurcating solution. Note that k=σ/δ≈π/2-σ/τ is close to π/2 and that κ=τ/δ is large since δ is small. Since k is close to π/2, 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 M is sufficiently big; analogous to the previous case we then obtain equation (3.32) and δ=σ/(π/2-σ/τ)2. For τ large, δ≈4σ/π2 and as before we have that k≈π/2-σ/τ is close to π/2 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 40% of the amplitude of the wave. We remark that the spatial shape of the approximation uϵ̂ given in (3.28) of the bifurcating solution to (2.1) is independent of the size of the bifurcation parameter μ>0 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 ϵ̂>0 go to zero, μ goes to zero and uϵ̂ approaches the null-solution—which is reasonable.Spreading of the wave. The automatic domain is located on [-1,1] and the FitzHugh-Nagumo parameters ϵ=0.05 and γ=0.5 are used. Depicted are the function ξ in (2.6) (giving the spatial shape of the approximative solution (3.28)) and the automaticity function in (2.2) for different values of α and β. The corresponding critical diffusion coefficient η0 and the values of k=(α/ϵ-γ)/η0 and κ=(β/ϵ+γ)/η0 in the function ξ are indicated as well. Upper row: increasing automaticity α of the automatic cells makes the wave spread further into the surrounding tissue. Lower row: decreasing stability β of the normal cells makes the wave spread further into the surrounding tissue.
(a)(b)(c)(d)(e)(f)


## 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 λ+≈iω0+λμμ, where λμ is the derivative of the principal eigenvalue with respect to μ, evaluated at μ=0. 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 μ>0. 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 η=η0-μ and calculate(3.2)0=ddμG=Gλλμ+Gμ
and then evaluate this expression at μ=0 to find λμ=-Gμ/Gλ. By using Maple an expression for λμ is found; denoting it by λμG it reads(3.3)λμG=κa3k2[(κ2+k2)sinh2(κa(M-1))+Mκ2](κ2+k2)sinh(2κa(M-1))+2aκ[(k2+κ2)cosh2(κa(M-1))-k2M](1-iγω0).Second, an alternative way to calculateλμ is given in [13]. Consider the spectral problem λζ=Fu(μ∣ζ), where at μ=0 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 μ=0 is solvable if and only if(3.4)λμ=〈Fuμ(0∣ϕ),ϕ*〉,〈ϕ,ϕ*〉=1.We turn to find the right-hand side〈Fuμ(0∣ϕ),ϕ*〉 in (3.4). The eigenfunctions ϕ and ϕ* occurring in (3.4) need to be normalized such that 〈ϕ,ϕ*〉=1. Let us start out by normalizing the eigenfunction ϕ given in (2.5) We could simply put c=1 in (2.5), but choose to determine it such that 〈ϕ,ϕ〉=1:(3.5)〈ϕ,ϕ〉=∫-MaMaϕ1ϕ¯1+ϕ2ϕ¯2dx=c2(1+ϵ)∫-MaMaξ2(x)dx,
and so the sought constant can be written as(3.6)c=1(1+ϵ)I0,
where we have defined(3.7)I0:=∫-MaMaξ2(x)dx.
Next, the complex constant d for the adjoint eigenfunction ϕ* in (2.17) is determined from the condition 〈ϕ,ϕ*〉=1. Using that ω0=1/ϵ-γ2 we calculate(3.8)〈ϕ,ϕ*〉=∫-MaMaϕ1ϕ¯1*+ϕ2ϕ¯2*dx=2cd¯ϵω0(ω0+iγ)I0
and find(3.9)d=12cI0(1+iγω0).
Now, Fuμ(0∣ϕ)=-(a2ϕ1,xx,0)=-(a2cξ′′,0) and(3.10)〈Fuμ(0∣ϕ),ϕ*〉=-∫-MaMaa2cξ′′(x)d¯ξ(x)dx=-cd¯a2I1=-a2I12I0(1-iγω0),
where I1 is the integral(3.11)I1:=∫-MaMaξ(x)ξ′′(x)dx.
So, the condition arising from (3.3) and (3.4) may be written as(3.12)λμG=-a2I12I0(1-iγω0).


## 3.2. Series Expansion

We attempt to find the solution to (2.10) as a power series in a small parameter ϵ̂:(3.13)(uϵ̂(x,s,ϵ̂)μ(ϵ̂)ω(ϵ̂)-ω0)=∑n=1∞ϵ̂nn!(un(x,s)μnωn).
Here uϵ̂ is the approximative solution, s=ω0t, μ is the bifurcation parameter, and ω is the time frequency.A Fredholm [17] solubility condition is given in [13], First, from 𝕁0u1=0 we may choose u1=z+z¯=eisϕ+e-isϕ¯. That is,(3.14)u1(x,s)=2cξ(x)(cos(s),ϵ(γcos(s)+ω0sin(s))).
It also holds that μ1=ω1=0 and(3.15)3(-iω2+μ2λμ)+3[Fuu(0∣u1∣u2),z*]+[Fuuu(0∣u1∣u1∣u1),z*]=0,
where z*=eisϕ*. However, in order to solve this equation for μ2 and ω2 we have first to solve(3.16)𝕁0u2=-Fuu(0∣u1∣u1),[u2,z*]=0
for u2, 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](3.17)iω2-μ2λμ=-2σ0σ1+4σ1σ¯1+2σ2σ¯2/3iω0,
where λμ is given by the right-hand side of (3.12) and(3.18)σ0=12〈Fuu(0∣ϕ∣ϕ),ϕ*〉,σ1=12〈Fuu(0∣ϕ∣ϕ¯),ϕ*〉,σ2=12〈Fuu(0∣ϕ¯∣ϕ¯),ϕ*〉.
We modify the equation to conform to the current series expansion and include the missing cubic term as well:(3.19)iω2-μ2λμ2=-2σ0σ1+4σ1σ¯1+2σ2σ¯2/3iω0+3ρ,
and here the cubic term reads(3.20)ρ=16〈Fuuu(0∣ϕ∣ϕ∣ϕ¯),ϕ*〉.
Recalling that Fuu(0∣ϕ∣ϕ)=2/ϵ((1-p(x))ϕ12,0) and Fuuu(0ϕϕϕ)=-6/ϵ(ϕ13,0), the coefficients in (3.19) may now be calculated. It is found that σ2=σ1=σ0 and(3.21)σ0=c2d¯ϵI2=I22I0ϵ(1+ϵ)I0(1-iγω0),ρ=-c3d¯ϵI3=-I32I02ϵ(1+ϵ)(1-iγω0),
where we have defined the integrals(3.22)I2:=∫-MaMa(1-p(x))ξ3(x)dx,I3:=∫-MaMaξ4(x)dx.
Now, solving (3.19) for μ2 and ω2 we obtain the bifurcation coefficient(3.23)μ2=22γI22-3ϵω02I0I3(ϵaω0I0)2(1+ϵ)I1.
The frequency coefficient ω2 is not significant for first-order accuracy but is given here for completeness:(3.24)ω2=-43I22(ω02+γ2)ω03(1+ϵ)I03ϵ2=-43I22(ω0I0ϵ)3(1+ϵ).


## 3.3. Summary and Discussion of Analytical Results

The bifurcation parameter for a small parameterϵ̂ is given by(3.25)μ=ϵ̂22γI22-3ϵω02I0I3(ϵaω0I0)2(1+ϵ)I1,
and the corresponding diffusion coefficient is δ=(η0-μ)a2 with η0 determined from (2.7). Here the frequency coefficient is ω0=1/ϵ-γ2 and the integrals are given by(3.26)I0=∫-MaMaξ2(x)dx,I1=∫-MaMaξ(x)ξ′′(x)dx,I2=∫-MaMa(1-p(x))ξ3(x)dx,I3=∫-MaMaξ4(x)dx
with p in (2.2) and the function ξ given in (2.6):(3.27)ξ(x)={coskx,|x|≤a,coskasinh(κa(M-1))sinh(κ(Ma-|x|)),a<|x|≤Ma
with k=(α/ϵ-γ)/δ0 and κ=(β/ϵ+γ)/δ0, δ0=η0a2. 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(3.28)uϵ̂(x,t)=2ϵ̂cξ(x)(cos(ω0t),ϵ(γcos(ω0t)+ω0sin(ω0t))),
where c=1/(1+ϵ)I0.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(3.29)λ+≈iω0+λμμ,λμ=-a2I12I0(1-iγω0),
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 a=1 and consider (2.7) with δ=η0. Now, put σ=α/ϵ-γ and τ=β/ϵ+γ such that (2.7) can be written as(3.30)σδ-arctan(τσ1tanh((M-1)τδ))=0.
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, M is big enough, such that we may approximate (3.30) with(3.31)σδ-arctan(τσ)=0.
Expanding the arctan-function yields the equation(3.32)σδ-(π2-στ)=0
on ignoring higher-order terms. Solving for the diffusion coefficient we get δ=σ/(π/2-σ/τ)2, which is small since σ is small. Consider now the function ξ above describing the spatial shape of the bifurcating solution. Note that k=σ/δ≈π/2-σ/τ is close to π/2 and that κ=τ/δ is large since δ is small. Since k is close to π/2, 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 M is sufficiently big; analogous to the previous case we then obtain equation (3.32) and δ=σ/(π/2-σ/τ)2. For τ large, δ≈4σ/π2 and as before we have that k≈π/2-σ/τ is close to π/2 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 40% of the amplitude of the wave. We remark that the spatial shape of the approximation uϵ̂ given in (3.28) of the bifurcating solution to (2.1) is independent of the size of the bifurcation parameter μ>0 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 ϵ̂>0 go to zero, μ goes to zero and uϵ̂ approaches the null-solution—which is reasonable.Spreading of the wave. The automatic domain is located on [-1,1] and the FitzHugh-Nagumo parameters ϵ=0.05 and γ=0.5 are used. Depicted are the function ξ in (2.6) (giving the spatial shape of the approximative solution (3.28)) and the automaticity function in (2.2) for different values of α and β. The corresponding critical diffusion coefficient η0 and the values of k=(α/ϵ-γ)/η0 and κ=(β/ϵ+γ)/η0 in the function ξ are indicated as well. Upper row: increasing automaticity α of the automatic cells makes the wave spread further into the surrounding tissue. Lower row: decreasing stability β of the normal cells makes the wave spread further into the surrounding tissue.
(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 α=0.03, β=0.1, ϵ=0.05, γ=0.5, and diffusion coefficient δ=(η0-μ)a2 is considered on [-Ma,Ma] with a=0.001 and M=2. From (2.7) we obtain that η0=0.053016 for purely imaginary eigenvalues ±iω0, ω0=1/ϵ-γ2=79/4.First, we check that condition (3.12) holds. The characteristic equation (2.8) enabled us to find an expression for λμG=-Gμ/Gλ. This expression is given in (3.3) and evaluating it for the current parameters gives λμG=0.82324-0.092621i. Denoting by λμ the right-hand side of (3.12) and checking the condition shows that |λμ-λμG|=0. Now, λ+=iω0+λμμ should agree to first order with λexact+ obtained from (2.8) with η=η0-μ. Note that since the real part of λμ is positive, the null-solution is unstable provided that μ>0. For different ϵ̂ we calculate μ given in (3.25) and check the errors |λ+-λexact+|. By using Müller's method, see [18], λexact+ is found to an accuracy of around 10-15. 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 [-Ma,Ma], M even (M=2 in the present case), with N+1 points, such that -a and a are centered between grid points. In this way the influence of the singularities at x=±a is diminished. It is readily found that we may choose N=M(2j+1), j=0,1,2,… to accomplish that ±a are located in the center of grid points. By discretizing (2.1) in space with homogeneous Dirichlet conditions we can formulate a semidiscrete system U′(t)=δDU(t)+H(U(t)). Here U=(V1,W1,…,VN-1,WN-1)T is the numerical solution vector and D is the 2(N-1)×2(N-1)-matrix of the classical second-order finite difference approximation of the second derivative acting on odd entries of U (that is on V). The function H reads(4.1)H2i-1(U)=1ϵf(Vi,Wi,pi),H2i(U)=Vi-γWi,
where f is given in (2.1) and i=1,…,N-1. We take N=1498 points and integrate the semidiscrete system in time up to T=4000 by the built-in  Matlab solver  ode15s with tolerance tol=10-12. 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:(4.2)u0=ϵ̂ccos2(πx2Ma)(1,ϵω02+γ2).We takeϵ̂=2·10-4 and obtain from Table 1 that μ=6.3081·10-4 and that the corresponding eigenvalue is λ+=5.2·10-4+4.44i. 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 c=28.833. Figure 2 depicts the approximative solution (3.28) and the numerical approximation in six snapshots covering 5/6th of the time-period 2π/ω0=1.4138. We observe a good agreement between the approximative solutions.Table 1
FHN eigenvalues. The rate corresponds to the convergence of λ+=iω0+λμμ towards λexact+ with respect to μ.
ϵ̂μλ+λexact+|λ+-λexact+|rate8.0000e-041.0093e-028.3089e-03 + 4.4432i8.4692e-03 + 4.4431i1.6235e-04—6.0000e-045.6773e-034.6738e-03 + 4.4436i4.7234e-03 + 4.4436i5.0241e-052.03854.0000e-042.5233e-032.0772e-03 + 4.4439i2.0869e-03 + 4.4439i9.7749e-062.01872.0000e-046.3081e-045.1931e-04 + 4.4440i5.1991e-04 + 4.4440i6.0555e-072.00641.0000e-041.5770e-041.2983e-04 + 4.4441i1.2986e-04 + 4.4441i3.7764e-082.00161.0000e-051.5770e-061.2983e-06 + 4.4441i1.2983e-06 + 4.4441i3.7738e-122.00021.0000e-061.5770e-081.2983e-08 + 4.4441i1.2983e-08 + 4.4441i3.6125e-162.0095Table 2
Convergence. The error is e=[uϵ̂-uref,uϵ̂-uref], where uϵ̂ is the analytical approximation and uref  is the numerical reference solution. The rate corresponds to the convergence with respect to ϵ̂.
ϵ̂μerrorerate8.0000e-041.0093e-028.0934e-05—6.0000e-045.6773e-034.8581e-051.77424.0000e-042.5233e-032.1892e-051.96592.0000e-046.3081e-045.2810e-062.0515FHN bifurcation solution. comparison of analytical and numerical approximations. Snapshots of the solutions at six different times; the analytical solution has been shifted in time to be in phase with the numerical solution. The dots indicate the automatic domain on [-0.001,0.001].
(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 Table1 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(4.3)[U,V]=12π∫02π〈U,V〉(s)ds,
where we recall that s=ω0t. 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.



---

*Source: 292183-2009-12-13.xml*