Stochastic dynamics of phase singularities under ventricular fibrillation in 2 D Beeler-Reuter model

The dynamics of ventricular fibrillation (VF) has been studied extensively, and the initiation mechanism of VF has been elucidated to some extent. However, the stochastic dynamical nature of sustained VF remains unclear so far due to the complexity of high dimensional chaos in a heterogeneous system. In this paper, various statistical mechanical properties of sustained VF are studied numerically in 2D Beeler-Reuter-Drouhard-Roberge (BRDR) model with normal and modified ionic current conductance. The nature of sustained VF is analyzed by measuring various fluctuations of spatial phase singularity (PS) such as velocity, lifetime, the rates of birth and death. It is found that the probability density function (pdf) for lifetime of PSs is independent of system size. It is also found that the hyper-Gamma distribution serves as a universal pdf for the counting number of PSs for various system sizes and various parameters of our model tissue under VF. Further, it is demonstrated that the nonlinear Langevin equat...


I. INTRODUCTION
Spatiotemporal phenomena in excitable media such as spiral waves and spiral wave turbulence have been studied extensively 1 in physics, chemistry, and biology.Studies on the spiral wave formation in cardiology are quite popular since the reentry of spiral wave plays an important role in the initiation of cardiac arrhythmias.Especially, spiral waves in the ventricular tissue may induce tachycardia, and the breakup of spiral waves leads to turbulent electrical activity called as ventricular fibrillation (VF).Elucidating the nature of VF in normal and diseased tissues may provide new insights into prevention and control of sudden cardiac death.
Many theoretical models have been developed to explain how spiral waves initiate and how a wave breaks in homogeneous and isotropic tissues with the use of simple qualitative cell models 2,3 and sophisticated ion channel cell models 4,5 (Ref.3 and references cited therein).Recent advances in computer technology have enabled us to perform simulations of three-dimensional, heterogeneous and anisotropic tissues with one of the sophisticated models. 6In addition, the initiation of VF in more realistic situations is now commonly investigated. 7On the other hand, the nature of VF in the steady-state has not yet been elucidated quantitatively.
VF is sustained by wandering reentrant waves wherein breaking and extinction of waves arise.These waves are classified into "spatiotemporal chaos" and "spiral wave turbulence": local shortlived excited reentrant waves could dramatically affect the evolving global spatiotemporal pattern.Spatial phase singularity (PS) 8,9 is an effective tool to quantify complex reentrant waves, which a akio.suzuki@aist.go.jp 2158-3226/2011/1(3)/032103/13 C Author(s) 2011 1, 032103-1 could be applied to data from experiments and numerical simulations.2][13][14][15] Gil et al. 11 investigated creation and annihilation mechanisms of defect pairs and predicted the probability distribution function of the number of defect pairs in oscillatory media.The prediction has been verified numerically and experimentally, 12 and modifications have been made to include the effect of boundary conditions 13 and noise. 14For spiral wave turbulence in 3D excitable media, different mechanisms are found and a new probability distribution function is proposed. 15For VF in the model tissues with many ion channels, the statistical properties of defects have not yet been elucidated in detail.
In this paper, we intend to build a stochastic dynamical model for describing the birth-death process of PSs in both normal and diseased tissues.For nonlinear complex phenomena like VF, one can introduce a stochastic model in the form of a nonlinear Langevin equation (NLE) or a corresponding Fokker-Plank equation. 16,17 here are two approaches to determine the form of NLE and to estimate the parameters of a model: One is based on a probability distribution, 18 and the other is based on time series data. 19ection II presents (A) Beeler-Reuter cardiac model in 2D tissue by using the method of numerical simulation; (B) a method of defining the phase with PS detection; (C) a hyper-Gamma distribution for the number of PSs, and a parameter estimation method for a corresponding nonlinear Langevin equation.In Section III, the statistical characteristics of PS under VF and the stochastic process of PS evolution are shown.In Section IV, detailed discussions are given on the statistical characteristics of VF in our model, system size dependency of VF, and validity of identified stochastic process.Section V is devoted to the conclusion.

A. Computational model
It is reported that electrical wave propagation in cardiac tissue can be described based on the following reaction diffusion equation: 20 where V m (mV) is transmembrane voltage, D (cm 2 /ms) is the diffusion coefficient, C m (μF/cm 2 ) is membrane capacitance, I ion (μA/cm 2 ) is the sum of ionic transmembrane currents, I st (μA/cm 2 ) is the stimulus current, and ∇ 2 is the Laplacian in 2D space.Various modifications can be made for describing ionic currents.Here, we have adopted the Drouhard-Roberge modification 21 of the Beeler-Reuter (BR) cardiac cell model 4 (BRDR model), which is the simplest ionic model for mammalian ventricular tissue using four ionic currents, five ionic gates, and intracellular calcium concentration.
Numerical simulations of the BRDR model are performed in a 2D square domain = [0, L] × [0, L] with Neumann (zero flux) boundary conditions.To reduce computation time, numerical integration is performed using a forward Euler scheme with a time step 0.025 ms.The Laplacian is discretized using a 9-point difference formula 20 with a spatial mesh 0.025 cm.Space-time variations of state variables are stored with a time step of 1.0 ms and a space step 0.05 cm.In our simulations, the values of the two parameters D = 9.72 × 10 −4 and C m = 1.0 are kept constant. 22The BRDR model parameter values of Ref. 4, 21 are adopted in this study, except for the values of slow inward current conductance g s (mS/cm 2 ).
Spiral waves are initiated by the same method as that described in the Ref. 20.First, a planar wave was initiated at the lower border by line stimulation over 2 ms at a strength of 40 μA/cm 2 , and it propagated upward.When a planar depolarization front reached to halfway across the 2D sheet, the state variables in the right half plane are set to the initial values.The wave front then curved and formed a spiral wave.4][25] In our numerical simulation, three slow inward current conductance values, g s = 0.09 (original value), 0.07 and 0.03 mS/cm 2 are examined.As g s decreases, the action potential duration (APD) of a cell also decreases as shown in Fig. 1(a).A similar nature can be seen for the wave forms in Fig. 1(b), wherein their wavelengths are 15.0, 11.9 and 4.1 cm, respectively.The dispersion relation (the relation between propagating velocity and stimulation frequency) is not changed and the slope of APD restitution curve flattens. 22or g s = 0.09 and 0.07 mS/cm 2 , an initiated spiral wave rotates a few times, and its breakup starts near the core 3 (Fig. 2(a)-2(f)), and then fibrillation develops.For g s = 0.03 mS/cm 2 , another type of breakup appears due to the Doppler effect of the meandering spiral wave 3 (Fig. 2(g)-2(i)).
The sustainability of fibrillation depends on the tissue size relative to the wavelength. 26If the tissue size is not large enough, fibrillation self-terminates quickly after the initiation of a spiral wave.Our simulation is therefore performed in extremely large domain (L = 25.0 cm).To assess the size effect, simulation is also performed in a relatively small domain (L = 20.0 cm for g s = 0.09 mS/cm 2 , L = 17.5 cm for g s = 0.07 mS/cm 2 and L = 12.5 cm for g s = 0.03 mS/cm 2 ).In each computation, the total duration of simulation is about 30.0 s from the planar wave initiation, and time series data over a stationary period (20.0 s for the period between 10.0 s and 30.0 s) are used in our statistical analysis.

B. Phase singularity
The state of VF may be characterized by a birth and death process of the phase singularity.To evaluate the number of phase singularities in a 2D system, one can use (i) the local finite-time Lyapunov exponent, 27,28 (ii) the Poincare-Hopf index, and (iii) a concise method based on a mutualphase between two state variables among the multiple variables in the BRDR model.Here we use the transmembrane voltage V m and the fast inward sodium current inactivation gate variable h.The local phase θ (x, y, t) at the point (x, y) is defined by  2) is shown as the angle from the horizontal vector and that one between the blue circle and the origin indicated by the red circle.
where arctan returns a value between −π and +π , (V * m , h * ) is set as the origin.To calculate this phase correctly, the origin must be encircled by all trajectories.In this study, the fixed origin (V * m , h * ) = (−50.0,0.01) is adopted for all simulations.These values are selected by observing the trajectories in the phase space.An example is depicted in Fig. 3.
The phase singularity (PS) is identified as the local site where the phase is undefined.The PSs appear at the core of the rotating waves, and their existence is a necessary (but not a sufficient) condition for the occurrence of breaking waves.The position of a PS can be identified from where θ is the local phase.The line integral is taken over r around a closed curve C surrounding a local site.The sign of Eq.(3) indicates the chirality; positive and negative signs correspond to clockwise and counterclockwise rotation, respectively.A spatial phase map was calculated for each frame of recorded data, and the phase singularities (PSs) are detected by checking the equality described above. 8The features of birth and death process of the PSs are analyzed (cf.Ref.10).A searching radius of 0.15 cm is adopted to link the PS from one frame to the next.The nearest PS in the next frame within the searching radius is linked.Unlinked PSs are regarded as ones having died or ones having just been born.Lifetime of a PS is estimated from information of its birth and death time.

C. Statistical model
As a statistical model of the PS evolution for the BRDR model, the Gamma and the hyper-Gamma (HG) distribution 18,29 are examined.The HG distribution (the probability density function (pdf)) is a generalization of the Gamma distribution for describing various complex phenomena in non-equilibrium open systems: where α and γ are the shape parameters and β is the scaling parameter.This distribution covers many different kinds of distributions (cf.Ref.18 in details).In the case of α = 1, for example, it reduces to the Gamma distribution.
There are infinite numbers of stochastic processes associated with the HG distribution.To model the statistical process of PSs, we choose the type II nonlinear Langevin equation (NLE) among the other types that give the same distribution as Eq.( 4) at the steady state: where m is a positive real number (m > 1) and the noise source η(t) is Gaussian-white η(t)η(t ) = 2Dδ(t − t ) with a null mean η(t) = 0.The parameters in Eq.( 4) are related to the ones in Eq.( 5) as All the parameters in Eq.( 6) (i.e., m, a, b and D) cannot be determined by estimating only the estimation of the parameters of the pdf (i.e., α, β, γ ).Thus we try to estimate the value of D directly from the time series data of PS.The general form of 1D nonlinear Langevin equation is written as where F(t) is Gaussian-white noise with the unit strength F(t)F(t ) = δ(t − t ) and a null mean F(t) = 0.The deterministic part g(x) and the stochastic part h(x, F(t)) are estimated by the relationships: 19 and Since the stochastic part becomes the value of noise strength D can be easily estimated.Then the values of the other parameters a, b, and m are estimated from the relations in Eq.( 6).

A. Statistical characteristics of PSs
Figure 4(a) shows the temporal evolutions of the number of PSs in a tissue with a domain size L = 25 cm for three values of slow inward current conductance g s = 0.09, 0.07, and 0.03 mS/cm 2 .Figure 4(b)-4(d) show the corresponding space-time plots of the PSs.The trajectories of the PSs are obviously different.In the case of g s = 0.09 and 0.07 mS/cm 2 , it seems that the PSs move rapidly and interact strongly with each other.In contrast, it seems that the movements of PSs are slow and their interactions are weak in the case of g s = 0.03 mS/cm 2 .To verify their significance in detail, the PS data in the VF state is statistically analyzed.The results are summarized in Table I.
First, let us look at the effect of g s in the L = 25 cm tissue.The mean number of PSs increases significantly, and the mean velocity (total traveling distance of a PS divided by its lifetime) decreases as the value of g s decreases.The notable points for the mean lifetime of PS and the mean birth and death ratio (the number of births and deaths per the number of PS) are as follows: (i) In the case of g s = 0.03 mS/cm 2 , the lifetime is five times longer, and (ii) the birth/death ratio is five times smaller than those of other cases.The nature of lifetime was investigated in detailed.The pdfs for lifetime are shown in Fig. 5. Exponential relaxation,P(t) ∝ exp(−at), appears in the long lifetime regions; a = 0.015 for g s = 0.09 mS/cm 2 , a = 0.018 for g s = 0.07 mS/cm 2 , and a = 0.002 for g s = 0.03 mS/cm 2 .Furthermore, several peaks appear in the short lifetime region for g s = 0.03 mS/cm 2 .Due to the significant difference, we can say that values of g s = 0.09 and 0.07 mS/cm 2 are large, whereas a value of g s = 0.03 mS/cm 2 is small.For larger values of g s , the mean lifetime, the standard deviation of PS density, and the number of deaths increase when the value of g s decreases.In contrast, the birth and death ratio decreases when the value of g s decreases.
Next, let us assess the size effect.The mean velocity, mean lifetime of PSs, mean birth ratio, mean death ratio in the inner region, and the pdf of lifetime are almost unchanged.These indices are size independence.The mean density of PS decreases and the standard deviation of the density of PS increases when the tissue size decreases.

B. Stochastic model of phase singularities
The suitability of stochastic model of PS evolution in fibrillated cardiac tissues is examined.The pdf (the Gamma and the hyper-Gamma distribution) of the number of PSs are calculated for all the cases.Their parameters are estimated by the maximum likelihood method.The theoretical probability densities are over-plotted with the ones from the numerical simulation of BRDR model as shown Fig. 6.The estimated parameters are listed in Table II.In the cases of large domain (L = 25 cm) simulations with g s = 0.07 and 0.03 mS/cm 2 , the Gamma distribution agrees well with the ones obtained from the simulations (Fig. 6(c) and 6(e)).In the other cases, however, the Gamma distributions do not agree with the ones obtained from simulations, especially in the tail regions (Fig. 6(a), 6(b), 6(d), and 6(f)).The obtained pdfs for the number of PSs show small probability in the right tail than that compared with the Gamma distribution.Consequently, it is found that the hyper-Gamma distribution can cover all the probability densities that are obtained in our simulation.For the system size L = 25 cm with g s = 0.07 and 0.03 mS/cm 2 , the identified hyper-Gamma distributions have profiles similar to the Gamma distribution.In the other cases, the identified hyper Gamma distribution describes the number of PSs well even in the tail region, where the Gamma distribution is not appropriate.Hence, it is quite appropriate to utilize the hyper-Gamma distribution to describe the number of PSs under VF in the BRDR tissues of different sizes with modified ion channel parameters.
Among the stochastic processes that give the hyper-Gamma distribution, the type II NLE in Eq.( 5) is considered to be the most suitable one.In order to confirm our prediction and to estimate the parameters, the relations in Eq.( 8) and Eq.( 9) are examined as shown in Fig. 7.It seems to be troublesome to estimate a deterministic part g(x) in Fig. 7(a) since data points are widely scattered.On the other hand, the linear relation, h(x) = a h x + b h , is observed in Fig. 7(b) in all cases.The linear coefficient a h is estimated by the standard least squares method.The parameter D of the   type II NLE is estimated from the relation in Eq. (10).The other parameters are determined using Eq.( 6).The estimated parameters are summarized in Table III.The noise intensity D becomes large when the tissue size decreases in all cases.
The Langevin simulations of the NLE with estimated parameters ( â, b, m and D) are carried out by the stochastic Runge-Kutta method 30 with a time step of 0.01 ms and an initial value x(0) = 1.0.Time evolutions for an L = 25 cm tissue is shown in Fig. 8.The means and the standard deviations from the Langevin simulation of the NLE with the estimated parameters listed Table IV can reproduce those obtained from the simulated PSs in the BRDR model.

IV. DISCUSSION
The statistical analysis of PSs, especially the probability density of the PS lifetime, indicates that at least two different types of VF can arise in our tissue model by modifying the slow inward current conductance g s in the BRDR cell model.Note here that the modification of the value of conductance may represent the normal and diseased tissues.The two different VFs sustained may have originated from different mechanisms of wave breakup.In the case of a large g s (normal or slightly damaged case), the reentrant wave breaks at a nearby PS and then the annihilation of two PSs, the original PS and one of the new-born PSs, occurs rapidly.Therefore the lifetime of PS has a tendency to have a short lifetime.On the other hand, in the case of small g s (heavily damaged case), the breakup of wave occurs at outside of the core.Hence the interaction between the original PS and new-born PSs is weaker and slower than the interactions among new-born PSs.A decrease in the mean birth/death ratio relative to the number of PSs supports our speculation.An important problem in the near future is how the exponential relaxation and observed peaks in the pdf of PS lifetime relate to the interactions of the PSs quantitatively.
The statistical analysis of PSs in tissues of several sizes shows that the mean of velocity, lifetime, birth ratio, death ratio in the inner region, and the pdf of the lifetime of PSs are independent of the tissue size and depend only on the cell parameter g s .We think that these are good indices to describe the time-averaged nature of stationary VF.On the other hand, the standard deviations of spatial PS  density and the pdfs of the number of PSs are strongly dependent on the tissue size.As the tissue size decreases, the standard deviations of spatial PS density increases in all cases.
For the pdf of the number of PSs, the relevant distribution changes from a Gamma-like distribution to a hyper-Gamma distribution.The case of g s = 0.09 mS/cm 2 is an exception.It seems that the tissue size L = 25 cm is not sufficiently large as compared to the propagating wave length.Although further investigation is required, one can say that the hyper-Gamma distribution gives a size-independent description for stationary VF.
Though there are many stochastic processes associated with the hyper-Gamma distribution, the type II NLE is selected as one of the appropriate processes.In the parameter estimation of the NLE, we neglect the constant term that appears in the estimated stochastic part of NLE.In spite of the simplification, the type II NLE can reproduce the temporal features of the PS for the BRDR model.
Gil et al. 11 have calculated the probability distribution of defects (phase singularities) based on the rate of creation and annihilation of defects in the framework of the Master equation approach.The effect of boundary condition, noise or the third dimensions are accounted in the subsequent papers [12][13][14][15] in the asymptotic regime.Their theoretical analyzes are based on the single state variable n (a pair of defects), the local nonlinear interactions involved in the system [11][12][13][14][15] can be categorized into the generalized Schlögl model: 31 X k 2 k 1 2X , 0 k 4 k 3 X and 2X k 6 k 5 0. In the present paper, we have studied the dynamics of defects (phase singularities: PSs) based on the rate of creation and annihilation of defects in the framework of the Langevin approach (or the Fokker-Planck approach).In this case, the features of both (a) local nonlinear interactions 31 with (b) global environmental fluctuation are accounted by introducing the fractional order nonlinearity and the multiplicative noise.Thereby, our model can mimic successfully both (i) the temporal order in the sense of time correlation function and (ii) the probability distribution of defects (PSs).
In our numerical simulations, the value g s of the conductance is modified in the BRDR model.Decreasing of inner ion currents has been reported for patients with congential long QT syndrome (LQT1 -LQT8). 32Especially, for the patients with LQT4 and LQT8 syndrome, the decreasing of I Ca , I Na and so on seems be relevant.It is also reported that I Ca increases in mild-to-moderated congestive heart failure (CHF), and it decreases in severe CHF. 33Acidosis and hypoxia in acute ischemic tissue tend to decrease of I Ca .On account of these facts, the modification of the g s for I Ca value is adopted as an initial trial of our study.

V. CONCLUSIONS
We have studied the stationary nature of VF in 2D BRDR tissues with normal and modified ionic current.It is shown that the time averaged characteristics of VF, such as velocity, lifetime, birth and death ratio, and the pdf of the lifetime of PSs are independent of the tissue size.It is also shown that the hyper-Gamma distribution can reproduce the probability distribution for the number of PSs shown for a wide variety of tissue size and modification of parameters (normal and diseased tissue).Further, the stochastic process associated with the type II nonlinear Langevin equation is shown to be an appropriate model for describing the number of phase singularities that arise during VF.The obtained statistical properties and the proposed stochastic process may provide new information for understanding of the nature of VF.

FIG. 2 .FIG. 3 .
FIG. 2. Breakup of a spiral wave in a 25 × 25 cm 2 tissue in BRDR model.Frames (a)-(c): normal current (g s = 0.09 mS/cm 2 ).The initiated spiral wave (a) starts to breakup close to the rotational tip (b) and VF is developed (c).Frames (d)-(f): modified current (g s = 0.07 mS/cm 2 ).Process is same as in normal tissue.Frames (g)-(i): modified current (g s = 0.03 mS/cm 2 ).The breakup of the spiral wave occurs at outside of the core due to the Doppler effect of meandering tip movement.

FIG. 6 .
FIG. 6. Probability density function of the number of PSs during VF in L×L cm 2 tissue with normal (g s = 0.09 mS/cm 2 ) and modified (g s = 0.07 and 0.03 mS/cm 2 ) slow inward current conductance.Blue lines indicate the theoretical probability density of Gamma distribution, and red lines indicate hyper Gamma distribution.(a) g s = 0.09 mS/cm 2 and L = 25.0 cm.(b) g s = 0.09 mS/cm 2 and L = 20.0 cm.(c) g s = 0.07 mS/cm 2 and L = 25.0 cm.(d) g s = 0.07 mS/cm 2 and L = 17.5 cm.(e) g s = 0.03 mS/cm 2 and L = 25.0 cm.(f) g s = 0.03 mS/cm 2 and L = 12.5 cm.

TABLE I .
Summary of statistical characteristics of phase singularities a number of birth or death per number of PS

TABLE II .
Estimated parameters of Gamma distribution P G (x; a, b) and hyper-Gamma P H G (x; α, β, γ ) distribution for the number of PSs

TABLE III .
Estimated parameters for type II nonlinear Langevin equation in Eq.(5) FIG.8.Time evolutions of the estimated type II nonlinear Langevin equations.Numerical integrations are carried out for 3 × 10 6 steps with a step size of 0.01, which corresponds to 30s for the VF simulations.The periods from 10 6 to 3 × 10 6 are shown.

TABLE IV .
The means and the standard deviations obtained from the nonlinear Langevin equation along with the estimated parameters.