Markov Chain Models of the Number of PS

0

…

Fig. 3.2 Illustration of the Markov process. For a given number of PS (NPS) the four possible transitions are shown. The only exception is for zero and one NPS where further removal of PS is either not possible at all or not possible through pair annihilation.

3.2 Defect Mediated Turbulence and Markov Chain Modelling

The properties of the cardiac dynamics during arrhythmia can be likened to those of other pattern forming systems such as hydrodynamic, turbulent ones. Also in turbu-lent systems topological defects occur which are comparable to the phase singularities (PS) observed in cardiac dynamics. Individual PS may behave similar to gas particles moving colliding and, unlike gas particles, randomly appearing. Viewing their paths and kinematics as, to a certain degree, random suggest the possibility of modelling the behaviour of the current number of PS (NPS) over time as a stochastic process.

For many turbulent systems the stochastic properties of the defect dynamics have been analyzed previously both experimentally and numerically [64–74]. The focus in these studies is to describe the stochastic properties of the number of defects and, using numerical models, to form a deeper understanding of the system properties.

Aside from the general descriptive nature of such a model, by comparison to numer-ical studies new insights may be gained from inspection of the model parameters. This includes insights into the underlying dynamical properties, which are difficult to quan-tify and understand in depth. This means that the approach may provide a method for better characterizing and learning about the dynamical properties of cardiac fib-rillation. After introducing the approach in the following sections, Sections 3.2.2 f.

provide a list of interpretations of the stochastic model, which may be applicable to the heart.

3.2.1 Markov Chain Models of the Number of PS

Two different approaches may be taken for characterizing the dynamics of PS, i.e.

either quantifying the statistics for only positive/negative PS individually or the sum of both positive and negative ones. Since typically the individual count is expected to be approximately half of the total count, a large difference between the two methods is not

3.2 Defect Mediated Turbulence and Markov Chain Modelling 27 expected. However, focusing on only a single chirality ignores the different mechanics of a PS entering compared to pair annihilation, and seems thus more accessible for simulations. This is because within a simulation framework with periodic boundary conditions, the topological conservation of charge means that PS always occur in pairs.¹

Here, I will derive the basic principles based on the assumption that the stochas-tic properties of the NPS can be described by a discrete, continuous time Markov jump model. In this context, the current state is described only by the probability distribution 𝑃 (NPS) and the basic transition can be split into four distinct rates or interactions which are illustrated in Figure 3.2:

𝐶(NPS) = 1

2𝑟_creation(NPS) 𝐸(NPS) = 𝑟_entering(NPS) 𝐴(NPS) = 1

2𝑟annihilation(NPS) 𝐿(NPS) = 𝑟_leaving(NPS).

(3.9)

The factor1/2 exists since each creation and annihilation event contributes two phase singularities to the rate of change. These rates thus describe the mechanisms:

𝐶(NPS): pair creation adds two PS.

𝐸(NPS): one PSentering the field of view or being created at a boundary.

𝐴(NPS): pair annihilation removes two PS.

𝐿(NPS): one PSleaving the field of view or colliding with a boundary.

These four different events are directly accessible from experiments when using track-ing methods. However, Qiao et al. [73] have suggested that an estimation from the NPS time series itself is possible even when the observation time step is large.²

One property of such a Markov model is the existance of a probability distribution 𝒫(NPS) which does not evolve over time. Given the above rates from Equation (3.9) conditions for the detailed balance in the case of stationary probability distribution 𝒫(NPS) can be derived to find it:

(𝐶(NPS) + 𝐸(NPS) + 𝐴(NPS) + 𝐿(NPS))⋅𝒫(NPS) =

𝐶(NPS− 2) ⋅ 𝒫(NPS− 2) +𝐸(NPS− 1) ⋅ 𝒫(NPS− 1) +𝐴(NPS+ 2) ⋅ 𝒫(NPS+ 2) +𝐿(NPS+ 1) ⋅ 𝒫(NPS+ 1)

(3.10)

1Since there are no boundaries with which PS can collide or at which they can be created.

2The idea here is, that the measured transition matrix is given by the multiplication of𝐷shorter time transition matrices: 𝒮(∆𝑡) = 𝒮(∆𝑡/𝐷) ⋅ 𝒮(∆𝑡/𝐷) ⋯. Thus for sufficiently large𝐷 (and thus terms), the actual rate can be estimated by solving the equation for𝒮(∆𝑡/𝐷)from the known𝒮(∆𝑡).

This calculating can be performed using the fractional matrix power.)

where𝒫(NPS)is set to zero for NPS< 0, since a negative NPS is not possible. Further, the transition probabilities can also be used to define the non-zero off diagonal elements of a transition matrix 𝒮. Finding the transition matrix provides a convenient way to describe the system of equations derived from Equation (3.10) and to solve it using symbolical or numerical tools.

For certain functional forms of the rates and limiting to a single family of positive or negative defects (NPS_±) analytical forms have been derived [65, 67]. Gil et al. [65]

describe that the annihilation rate in such systems is often quadratic in NPS_± leading to a squared Poisson distribution for𝒫(NPS_±). Further, Daniels et al. [67] derive the stationary distribution for a linear leaving rate, a quadratic annihilation rate, and a constant creation and entering rate.

A general quadratic model with the constants const, 𝑚, and𝑎:

𝑟_{all rates}(NPS) =const+ 𝑚 ⋅NPS+ 𝑎 ⋅NPS² (3.11) will be used in many parts for all rates in this thesis. Although limitations to the model will be used when less data is available. In this case the rates are given by:

𝑟_{creation,entering}(NPS) =const+ 𝑚 ⋅NPS (3.12) 𝑟annihilation,leaving(NPS) = 𝑚 ⋅NPS+ 𝑎 ⋅NPS² (3.13) These limitations disregard the quadratic term for the entering and creation rates, and the constant term for the leaving and annihilation rate. In both cases these terms are theoretically not expected to be important or, in the latter case, expected to be zero³.

In the results Chapter 6 the stationary distribution𝒫(NPS)for the rates described by Equation 3.11 will be found by symbolically or numerically solving Equation (3.10) with a sufficiently high maximum NPS. Further properties of the Markov model system can be readily estimated numerically by simulating the system using the Gillespie algorithm [75].

Alternative Models for the Stationary Distribution

As an alternative to the Markov model derived above the following models will be considered for the description of the stationary probability distribution of the NPS:

• As a generic model and approximation expected to be correct for large enough NPS a Gaussian model can be fitted to the data. The Gaussian model has been used by Sugimura and Kori [76] for numerical simulations to estimate the time until the chaotic dynamics stop – a property commonly seen in such systems [77]. Additionally to a Gaussian probability distribution the authors introduce a transition rate from the state with one PS to the actual quiescent

3The annihilation/leaving rate has to be zero for NPS = 0. Any negative rate arising due to fitting or extrapolation is thus replaced with zero.

3.2 Defect Mediated Turbulence and Markov Chain Modelling 29

Fig. 3.3 Entering (𝐸), leaving (𝐿), creation (𝐶) and annihilation (𝐴) rates as function of the positive defects in an inclined-layer (30°) convection as presented by Huepe et al. in [70].

𝜖is defined by the temperature ratio 𝜖 = Δ𝑇 /Δ𝑇_𝑐− 1and is a parameter tuning the chaos in the system. Reprinted from [70], with the permission of AIP Publishing.

state where no activity is left. The approach of using a Gaussian distribution is motivated by the central limit theorem.

• A Markov model that only allows for a single step increment or decrement in the NPS. A comparison with such a model highlights the importance of pair interactions for the dynamics.