A stochastic phase model with reflective boundary and induced beating for the cardiac muscle cells
aa r X i v : . [ n li n . AO ] J u l A STOCHASTIC PHASE MODEL WITH REFLECTIVE BOUNDARY ANDINDUCED BEATING FOR THE CARDIAC MUSCLE CELLS
GUANYU ZHOU , TATSUYA HAYASHI AND TETSUJI TOKIHIRO Abstract.
We consider the stochastic phase models for the community effect of cardiac muscle cells. Themodel is the extension of the stochastic integrate-and-fire model in which we incorporate the irreversibilityafter beating, induced beating and refractory. We focus on investigating the expectation and varianceof (synchronized) beating interval. In particular, for the single-isolated cell, we obtain the closed-formexpectation and variance of the beating interval, and we discover that the coefficient of variance (CV)has upper limit p /
3. For two-coupled cells, we derive the partial differential equations (PDEs) for theexpected synchronized beating intervals and the distribution density of phase. Moreover, we also considerthe conventional Kuramoto model for both two- and N -cells models, where we establish a new analysis usingstochastic calculus to obtain the CV of the “synchronized” beating interval, and make some improvementto the literature work [21]. Introduction
A cardiac muscle cell (cardiomyocyte) has a distinguishing property among biological cells; it generatesspontaneous pulsation. Heartbeat is a macroscopic phenomenon in which pulsations of cardiac musclecells are tuned to a certain rate. Since each cell has its own beating rhythm when isolated, there mustbe a certain mechanism for synchronizing pulsations of cardiac muscle cells. Extensive works have beendevoted to understanding this mechanism both experimentally and theoretically [1, 6, 9, 10, 11, 16, 24,26, 27, 29, 33, 37, 39]. Contraction of a cardiac muscle cell is caused by complex electrophysiologicalprocesses and detailed analyses require elaborated mathematical models composed of a huge number ofequations [13]. To understand the essence of synchronization, however, a small number of simultaneousordinary equations of membrane currents and action potentials, such as the Hodgkin-Huxley equation[15] or its reduced forms, the FitzHugh-Nagumo (FN) equation [8, 31] and the Van der Pol equation (cf.[18, 30]), are enough to capture the key phenomenon of the cell dynamics.The cardiac muscle cells in a tissue are individual entities with identical genetic informations; however,these difference of individual cells are ironed out when they becomes clusters or tissues, which is calledthe “community effect” of cells as induced uniformity [17, 20]. Besides of the individual information(for example, the dynamics of the membrane currents of individual cells), to achieve a comprehensiveunderstanding of the cardiomyocytes’ dynamics, the analysis of the epigenetic information (the communityeffect) is mandatory. Since it is difficult to control the conditions and qualities of cells, there existslimitations in the biological experiments to study the community effect. To overcome these problem,the mathematical modeling is one of the most powerful approaches. In the present paper, we intend tounderstand the community effect of cardiomyocytes by proposing and studying the mathematical models,which should incorporate the essential properties of the biological system, and can somehow reproducethe experimental results [17, 20].To investigate the community effect of cardiomyocytes, we modify the conventional Kuramoto model[21, 22] by incorporating the conceptions of irreversibility of beating, the induced beating and refractoryto capture the essential properties of cardiomyocytes’ synchronization. Our model can be regarded as anmodification of the stochastic phase model or the integrate-and-fire model [4, 28], which has been widelyused as a spiking neuron model [2, 19, 32, 34]. We utilize the phase models for two reasons. First, fromthe biological experiments [17, 20], only the data of beating intervals is available. However, the Hodgkin-Huxley, FitzHugh-Nagumo or Van der Pol equations [15, 8, 31, 18, 30] model the dynamics of membranecurrents or ion concentration. Without adequate information of potential and ion concentration, it ishard to determine the parameters of these equations appropriately. Moreover, the application of these , TATSUYA HAYASHI AND TETSUJI TOKIHIRO models to each cell in N -cells network ( N ≫
1) yields a large number of nonlinear equations, whichis difficult to dealt with. Next, since we mainly focus on investigating the beating intervals, one canthink of the cardiomyocyte with rhythmic beating as an oscillator. Then, the phase models [21, 22, 38]are suitable mathematical tools to analyze the oscillation. In fact, the stochastic phase model is wellapplicable to model the distribution of the beating intervals (oscillation periods). For instance, in thecase of single-isolated cell (single oscillator), one can decide the intrinsic frequency and noise strength ofthe phase model by the experimental data of beating intervals and the formulas (3.3), (3.5). In addition,one can derive the phase equation from the FN model (see [22]).The main contribution of this paper is summarized in two aspects. First, it is an original idea toincorporate the stochastic phase equation with reflective boundary, induced pulsation and refractory, tomodel the (synchronized) oscillation of cardiomyocytes. [14] compares the simulation of the proposedmodels with the observation from biological experiments [17, 20], which indicates the well applicabilityof our models. The present paper, as a theoretical supplement to [14], is only devoted to the theoreticalanalysis. For single-isolated cardiomyocyte, we obtain the explicit relationship between the parameters(intrinsic frequency and noise strength) of the model and the statistic properties (expectation and vari-ance) of the beating interval. For two-coupled cells, by the renewal theory and Fokker-Planck equation,we derive the PDEs associated with the expectation the synchronized beating intervals and the distribu-tion density of phases. Although we cannot obtain the closed-form of the statistic properties, the PDEswith non-standard boundary conditions deserve the comprehensive theoretical/numerical analysis fromthe mathematical points of view.Second, we also consider the conventional phase model, and make several improvement to the existingresults [21]. In particular, we present a rigorous calculation of the coefficient of variance (CV) for bothtwo- and N -cells models using the theories of Itˆo integral, thanks to which, we provide the formulas todetermine the proper reaction coefficients of the model for the case of two-coupled cells.The rest of this paper is organized as follows. In Section 2, we introduce the biological background ofour work, and explain the connection between the FN model and the phase eqaution. Section 3 is devotedto the model with reflective boundary for the single-isolated cell. We study the stochastic phase modelsfor two-coupled cells in Section 4. The N-cell network is dealt with in Section 5. The concluding remarkis addressed in Section 6.2. From the FitzHugh-Nagumo model to the phase model
As a preliminary, we briefly introduce the The experimental approach to understand the epigeneticinformation of cardiomyocytes. And then, let us explain the connection between the the FN model andthe phase model.2.1.
The experimental approach.
The on-chip cellomics technology has been applied to investigatingthe community effect of cardiomyocytes [17, 20], which, simply speaking, includes three steps: (1) Thecells are taken from a community/tissue using a nondestructive cell sorting procedure. (2) We put the cellsin a microchamber (on chip) where we can design the cell network and control the medium environment.(3) We measure the beating intervals of each cell on chip by light signal (not the membrane currents).The procedure of the bio-experiment is described in Figure 2.1 (a) (see [17, 20] for details).To analyze the distribution of the beating intervals from experiments (see Figure 2.1 (right)), one canapply the FN equations with noise (2.12) to model the dynamics of the membrane currents. For example,in Figure 2.3 (c)(d), we plot the trajectories of FN model with rhythmic action potential influenced bynoise. However, it is nontrivial to determine the suitable parameters for FN model such that the distri-bution of beating interval generated by simulation (see Figure 2.4 (b)) coincides with the experimentaldata (Figure 2.1 (b)) well. To tackle this problem, we regard the cardiomyocyte as a oscillator satisfyingthe phase model (2.13) with intrinsic frequency µ and noise strength σ . This two parameters ( µ, σ ) canbe easily determined from the bio-experimental data of the beating interval. For above reason, we utilizethe phase model instead of FN model. In fact, the phase model can be derived from the FN system.2.2. From the FitzHugh-Nagumo model to the phase model.
The FN model has been widelyapplied to model the membrane current of the spiking neuron or cardiomyocyte, which can be regarded
STOCHASTIC PHASE MODEL FOR CARDIAC MUSCLE CELLS 3 (a) (b)
Figure 2.1. (a) The on-chip cellomics technology. (b) An example of the experimentaldata of the beating interval of two cardiomyocytes before and after coupling.as a simplification of the famous Hodgkin-Huxley model. First, let us pay attention to the case of thesingle cell without noise effect, the FN model of which is given by: dudt = u ( u − a )(1 − u ) − w, (2.1a) dwdt = τ ( u − bw ) . (2.1b)Here u denotes the membrane current, and ( a, b, τ ) the parameters ( τ ≪ w decreases below 0, u increase instantly, which corresponds to the pulsation of the membrane potential (beating). Therefore, onecan regard that w is associated to the refractory ( w also depends on u ). For a = − . b = 0 . τ = 0 . u behaves like a T -periodic function (see Figure 2.2 (a) with T ≈ t , the trajectory ( u, w ) tends to a limit cycle , that is, ( u ( t + T ) , w ( t + T )) =( u ( t ) , w ( t )) (see Figure 2.2 (b)). Hence, one can find a homeomorphism which maps the points ( u ( t ) , w ( t ))on the limit cycle to the phase function φ ( t ) given by(2.2) dφ ( t ) = µdt, where µ = T denotes the intrinsic frequency. Here, φ is also T -periodic if we set φ = φ + k ( ∀ k ∈ Z ), orequivalently φ takes value in torus [0 ,
1) (see Figure 2.2 (d)), which means φ jumps to 0 when approaching1 ( φ ( t ) = 0 when φ ( t − ) := lim s ↑ t φ ( s ) = 1). See Figure 2.2 (c) for an example of φ .Denote by x ( t ) := ( u ( t ) , w ( t )) the trajectories of system (2.1), and by χ ( t ) := ( u ( t ) , w ( t )) the trajec-tories of limit cycle, i.e., χ ( t + T ) = χ ( t ). Then, since φ is T -periodic, we can regard φ ( t ) as a functionof χ , i.e., φ ( t ) = φ ( χ ( t )). In fact one can extend such φ to all trajectories x , namely φ ( x ( t )) (see [22] forthe detailed argument).Setting f ( x ) := [ u ( u − a )(1 − u ) − w, τ ( u − bw )] ⊤ , we find that(2.3) dφ ( x ( t )) dt = ∂φ∂ x · d x dt = ∂φ∂ x · f ( x ) . Putting together with (2.2),(2.4) µ = dφ ( x ( t )) dt (cid:12)(cid:12)(cid:12)(cid:12) x = χ ( φ ) = ∂φ∂ x (cid:12)(cid:12)(cid:12)(cid:12) x = χ ( φ ) · f ( χ ( φ )) . In brief, to model the dynamics of the membrane currents, one can apply the FN model involving withtwo variables ( u, w ) and parameters ( a, b, τ ). Meanwhile, to describe the rhythmic oscillation, the phasemodel with intrinsic frequency µ (or the beating interval T ) is sufficient. GUANYU ZHOU , TATSUYA HAYASHI AND TETSUJI TOKIHIRO (a) (b) (c) (d) Figure 2.2. (a) The periodic solution u ( t ) of the FitzHugh-Nagumo model (2.1). (b) Theperiodic trajectories of ( u ( t ) , w ( t )). (c)(d) The phase model (2.2) with φ in torus [0 , u and u represent the membrane currents for two cells respectively, whichsatisfies the coupled FN model: du dt = u ( u − a )(1 − u ) − w + κ ( u − u ) , dw dt = τ ( u − bw ) , (2.5a) du dt = u ( u − a )(1 − u ) − w + κ ( u − u ) , dw dt = τ ( u − bw ) . (2.5b)Here, κ ( u i − u j ) describes the interaction between two cells. In Figure 2.3 (a)(b), we show an example ofthe synchronization of ( u , u ), where the trajectories ( u , w ) and ( u , w ) tend to the same limit cyclefor sufficiently large time t .Assume that trajectories { ( u i , w i ) } i =1 , are synchronized and T -periodic for large t , that is, { ( u i , w i ) } i =1 , both tend to the limit cycle χ with χ ( t + T ) = χ ( t ) ( i = 1 , x i ( t ) := ( u i ( t ) , w i ( t )) , κ ( x j − x j ) := [ u j − u i , ⊤ ( i, j = 1 , , i = j ) , and rewrite the model (2.5) as follows d x dt = f ( x ) + κ ( x − x ) , (2.6a) d x dt = f ( x ) + κ ( x − x ) . (2.6b)Denoting by φ ( t ) and φ ( t ) the phase functions for x and x respectively, we consider φ i ( t ) as afunction of the limit cycle χ , i.e., φ i ( t ) = φ i ( χ i ( t )), satisfying(2.7) dφ i dt = µ ( µ = 1 T ) . Here, φ i is T -periodic if we set φ = φ + k ( ∀ k ∈ Z ). Reversely, one can think of χ as a function of φ i , saying χ ( φ ( t )). Note that φ i ( χ ( t )) can be extended to all trajectories x i , namely φ i ( t ) = φ i ( x i ( t )).Analogously to (2.3) and (2.4), on the limit cycle χ ( φ i ),(2.8) µ = dφ i ( x i ( t )) dt (cid:12)(cid:12)(cid:12)(cid:12) x i = χ ( φ i ) = ∂φ i ∂ x i (cid:12)(cid:12)(cid:12)(cid:12) x i = χ ( φ i ) | {z } =: Z ( φ i ) · [ f ( χ ( φ i )) + κ ( χ ( φ ) − χ ( φ ))] = Z ( φ i ) · f ( χ ( φ i )) , where κ ( χ ( φ ) − χ ( φ )) = κ ( χ ( µt ) − χ ( µt )) = 0 because the synchronization φ = φ = µt occurs at χ .And Z ( φ i ( χ ( t ))) is also T -periodic.Under the assumption that the difference between x i and χ are small, that is | x i − χ | = O ( ε ) ≪
1, wecalculate as(2.9) dφ dt = ∂φ ∂ x · [ f ( x ) + κ ( x − x )]= ( Z ( φ ( t )) + O ( ε )) · [ f ( χ ( φ )) + κ ( χ ( φ ) − χ ( φ )) + O ( ε )]= µ + Z ( φ ( t )) · κ ( χ ( φ ) − χ ( φ )) + O ( ε ) . STOCHASTIC PHASE MODEL FOR CARDIAC MUSCLE CELLS 5 (a) (b) (c) (d)
Figure 2.3. (a) The synchronization of u ( t ) and u for the two-coupled FitzHugh-Nagumo models (2.5). (b) The periodic trajectories { ( u i ( t ) , w i ( t )) } i =1 , have the samelimit cycles. (c) The periodic solution u ( t ) of the FitzHugh-Nagumo model for one cellwith noise (2.13). (d) The periodic trajectories ( u ( t ) , w ( t )) of (2.13).Let us adopt the approximation approach from [22]. For sufficiently large time t , ( u , w ) and ( u , w )tend to the limit cycle. Then Z ( φ ( t )) · κ ( χ ( φ ) − χ ( φ )) can be regarded as a perturbation term, whichis approximately replaced by its average in ( t, t + T ), that is1 T Z t + Tt Z ( φ ( t ′ )) · κ ( χ ( φ ( t ′ )) − χ ( φ ( t ′ ))) dt ′ = 1 T Z T Z ( φ ( t ) + µt ′ ) · κ ( χ ( φ ( t ) + µt ′ ) − χ ( φ ( t ) + µt ′ )) dt ′ + O ( ε )= Z T Z ( η + φ ( t ) − φ ( t )) · κ ( χ ( η ) − χ ( η + φ ( t ) − φ ( t ))) dη =: Γ( φ ( t ) − φ ( t )) , where we have ignored the small term O ( ε ). Replacing Z ( φ ( t )) · κ ( χ ( φ ) − χ ( φ )) by the averageΓ( φ ( t ) − φ ( t )) in (2.9) yields(2.10) dφ dt = µ + Γ( φ ( t ) − φ ( t )) . We have presented a rough derivation of the phase equation for φ above ( φ can be treated in the sameway). The obtention of the closed-form of Γ( · ) requires technical calculation, which is omitted here. Onecan refer to [22] for more rigorous and detailed mathematical arguments. For simplicity, we replace Γ( · )by sin(2 π ( · )) without losing the essentiality of the model.In summary, we chose the Kuramoto model as the basic model to study the synchronization behaviorof two-coupled cardiomyocytes:(2.11) dφ dt = µ + sin(2 π ( φ ( t ) − φ ( t ))) ,dφ dt = µ + sin(2 π ( φ ( t ) − φ ( t ))) . Remark 2.1.
We apply the Kuramoto model owing to its wide application in studying the oscillators’synchronization. But from the mathematical points of view, it is worth to consider other interaction termsbesides of sin(2 π ( φ j − φ i )), for example ( φ j − φ i ) / | φ j − φ i | , ( φ j − φ i ) α and so on.2.3. The FitzHugh-Nagumo and phase models with noise.
The biological experiments (Figure 2.1(b)) show that the beating intervals of cardiomyocytes are not perfectly periodic, which indeed are effectedby noise. Therefore, it is necessary to consider the FN model with noise: du = [ u ( u − a )(1 − u ) − w ] dt + σ dW ( t ) , (2.12a) dw = τ ( u − bw ) dt, (2.12b) GUANYU ZHOU , TATSUYA HAYASHI AND TETSUJI TOKIHIRO (a) (b) (c) Figure 2.4. (a) A realization of the phase φ from (2.12). (b) The distribution ofthe beating intervals ∆ t from the FitzHugh-Nagumo model (2.12) with ( a, b, τ, σ ) =( − . , . , . , . µ = E (∆ t ) and σ = p Var (∆ t ) µ according to the formula (2.14).where σ > dW ( t ) the white noise ( W ( t ) is the standard Brownianmotion). A realization of the membrane current u ( t ) and the trajectories ( u ( t ) , w ( t )) of (2.12) is presentedin Figure 2.3 (c)(d). A simulation of the distribution of the beating interval ∆ t is plotted in Figure 2.4(b)) with mean value E (∆ t ) ≈ .
88 and standard variance p Var (∆ t ) ≈ . a, b, τ, σ ) ofFN model has not been understood fully, the phase model with noise is applicable to study the beatingprocess of cardiomyocyte, which is stated as follows: dφ ( t ) = µdt + σdW ( t ) , (2.13a) φ (0) = 0 , (2.13b)where σ > t := t ( k ) − t ( k − with t ( k ) the first passage time that φ ( t ) = k . Orequivalently, we set φ ( t ) = 0 when φ ( t − ) = 1 (the phase jumps to 0 when reaching 1), and ∆ t = t ( k ) − t ( k − with t ( k ) := inf { t > t ( k − | φ ( t − ) = 1 } ( t (0) = 0). By stochastic calculus, one can verify that(2.14) E (∆ t ) = 1 µ , Var (∆ t ) = σ µ , CV (∆ t ) = s Var (∆ t ) E (∆ t ) = σ √ µ . By (2.14), together with E (∆ t ) = 108 .
88 and p Var (∆ t ) = 17 .
662 (the simulation from FN model),we first compute the parameter ( µ, σ ) = (0 . , . u ( t ) in Figure 2.3 (c), when u ( t ) increases from 0 . . − . u , and also tothe period of oscillation cycle. In other words, when the action potential (the pulsation of u , or calledbeating) occurs, the noise effect is somehow inhibited such that the pulsation cannot be reversed by thenoise. This irreversibility has not been captured by the phase model (2.13), which may be the main reasoncausing the inconsistency between the distributions of Figure 2.4 (b) and (c). To address issue, in thenext section, we will propose a phase model incorporating the irreversibility after beating.3. The phase models for an isolated cell
Regarding the single-isolated cardiomyocyte as an oscillator with phase φ , we say the cell beats at time t when φ ( t − ) = 1, and we let φ ( t ) jumps to 0 immediately after beating to start a new oscillation cycle(beating interval). In view of the dynamics of u in Figure 2.3 (c), the beating ( φ ( t − ) = 1 , φ ( t ) = 0) STOCHASTIC PHASE MODEL FOR CARDIAC MUSCLE CELLS 7 corresponds to the action potential (i.e., u increases quickly from 0 . φ ( t ) = 0 such that φ cannot bedragged backward by noise. To incorporate the irreversibility after beating, we consider the stochasticphase model with the reflective boundary : dφ ( t ) = µdt + σdW ( t ) + dL ( t ) , (3.1a) φ (0) = 0 , (3.1b) φ ( t ) = 0 when φ ( t − ) = 1 , (3.1c)where L ( t ) is the process to prevent φ ( t ) being driven backward by noise when φ ( t ) = 0. L ( t ) indeeddescribes the reflective boundary at φ = 0 (cf. [12, 23, 35]). To state the definition of L ( t ), we set t ( k ) the k -th time that φ ( t ) approaches 1 ( k = 1 , , . . . ):(3.2) t ( k ) = inf { t > t ( k − : φ ( t − ) = 1 } , t (0) = 0 . We say the cell beats at t ( k ) ( k = 1 , , . . . ). In view of (3.1c), when φ approaches 1 (i.e. φ ( t ( k ) − ) = 1), φ immediately jumps to 0, i.e., φ ( t ( k ) ) = 0, and a new oscillation cycle begins.The rigorous definition of L ( t ) is given as follows (see Figure 3.1 (a)):(L1) L (0) = 0, and L ( t ) is a nondecreasing, continuous process for t ( k − ≤ t < t ( k ) such that φ ( t ) ≥ L ( t ) increases only when φ ( t ) = 0.In simulation, a simple approach [36] is to reset φ ( t n ) = 0 when φ ( t n ) = φ ( t n − ) + µ ∆ t + √ ∆ tση < t denotes the time-step (∆ t ≪ t n = n ∆ t ( n = 1 , , · · · , ) and η ∼ N (0 , Remark 3.1.
Assume that the phase state φ (0 < φ <
1) corresponds to the promptly decreasing ofthe membrane current u from 0 . − .
1, which also cannot be reversed by noise. And one may alsoinput a reflective boundary at φ = φ to control the noise, such that φ cannot be driven backward bynoise when φ = φ . But in this paper, we only implement the reflective boundary at φ = 0, which seemssufficient for application [14].During the first beating process, i.e., t ∈ [0 , t (1) ), integrating (3.1a) yields φ ( t ) = µt + σW ( t ) | {z } =: X ( t ) + L ( t ) , which means φ ( t ) is in fact a combination of the ( µ, σ )-Brownian motion X ( t ) and the process L ( t ). Weplot an example of ( φ ( t ) , X ( t ) , L ( t )) for t ∈ [0 , t (1) ) in Figure 3.1 (a). Since φ returns to 0 instantly whenapproaching 1, φ ( t ) is a renewal process for t ≥ t ( k ) ( k = 1 , , · · · ) (see Figure 3.1 (b)), and the beatingintervals (oscillation period) ∆ t ( k ) := t ( k ) − t ( k − ( k = 1 , , . . . ) are independent, identically distributedrandom variables. Hence, we only need to investigate ∆ t (1) = t (1) − Theorem 3.1.
For µ ≥ , σ > , we have (3.3) E ( t (1) ) = ( /σ for µ = 0 , ( θ − e − θ ) / ( µθ ) for µ > , (3.4) Var ( t (1) ) = ( / (3 σ ) for µ = 0 ,e − θ ( e − θ − e θ + 2 θe θ + 4 + 4 θ ) / ( µ θ ) for µ > , where θ = 2 µ/σ . Remark 3.2.
Throughout this paper, we always consider the non-negative intrinsic frequency ( µ ≥ σ > µ →
0, one can validate θ − e − θ µθ → σ and e − θ ( e − θ − e θ +2 θe θ +4+4 θ ) µ θ → σ . The coefficient of variance (CV) is given by:(3.5) CV ( t (1) ) = Var (∆ t ) E (∆ t ) = ( / µ = 0 ,K ( θ ) for µ > , GUANYU ZHOU , TATSUYA HAYASHI AND TETSUJI TOKIHIRO (a) (b) (c) Figure 3.1. (a) A realization of ( φ ( t ) , X ( t ) , L ( t )) with ( µ, σ ) = (0 . , X ( t ) is the( µ, σ )-Brownian motion, φ ( t ) = X ( t ) + L ( t ) and L ( t ) increases only when φ ( t ) = 0 suchthat φ ≥ φ ( t ) with ( µ, σ ) = (0 . ,
1) by numericalsimulation [36]. (c) The distribution of the beating interval t (1) of (3.1) with the identicalmean value and variance of the FitzHugh-Nagumo model.where K ( θ ) = e θ ( e − θ − e θ +2 θe θ +4+4 θ )(1+ θe θ − e θ ) and K ( θ ) ↑ / θ ↓ Remark 3.3.
We can compute the mean value and CV of the beating interval from the experimental data.From (3.5), we first determine θ . Then, with θ and the mean value, the coefficient ( µ, σ ) can be calculatedby (3.3). The applicability of our model (3.1) to the bio-experimental data has been discussed in [14].Here, we only compare the simulation with the FN model presented in Section 2. Using E (∆ t ) = 108 . p Var (∆ t ) = 17 .
662 from the simulation of the FN model (see Figure 2.4 (b)), we calculate thecorresponding ( µ, σ ) of (3.1), and carry out the simulation to plot the distribution of the beating interval(see Figure 3.1 (c)).
Remark 3.4.
Besides of the reflective boundary, the phase-dependent noise strength σ ( φ ) with σ (0) = 0is also considerable. In fact, this approach has been applied to the Langevin equation with noise modelingthe ion channels’ dynamic [5] for the Hodgkin-Huxley formulation, where the noise effect is inhibited whenthe proportion of the opened channels reaches 0 or 1 such that the proportion is bounded in [0 , Proof of Theorem 3.1.
For any function g ( x ) in [0 ,
1] with continuous differential ∂g∂x and ∂ g∂x , Ito’s formulayields (0 ≤ t ≤ t (1) ):(3.6) g ( φ ( t − )) = g ( φ (0)) + Z t (cid:20) µ ∂∂x + σ ∂ ∂x (cid:21) g ( φ ( s )) ds + Z t σ ∂g∂x ( φ ( s )) dW ( s ) + Z t ∂g∂x ( φ ( s )) dL ( s ) . Since L ( t ) is nondecreasing and increases only when φ ( t ) = 0 (see (L1)(L2)),(3.7) Z t ∂g∂x ( φ ( s )) dL ( s ) = Z { ≤ s ≤ t | φ ( s )=0 } ∂g∂x (0) dL ( s ) . Now, let g ( x ) be the solution of: (cid:20) µ ∂∂x + σ ∂ ∂x (cid:21) g ( x ) = − < x < , (3.8a) g (1) = 0 , ∂g∂x (0) = 0 . (3.8b) STOCHASTIC PHASE MODEL FOR CARDIAC MUSCLE CELLS 9
In view of φ (0) = 0 and φ ( t (1) − ) = 1, from (3.8), (3.7) and (3.6), we obtain(3.9) 0 = g (0) + Z t (1) − ds | {z } = − t (1) + Z t (1) σ ∂g ∂x ( φ ( s )) dW ( s ) . Because the expectation of an Itˆo’s integral is zero ([7, 25]), E "Z t (1) σ ∂g ∂x ( φ ( s )) dW ( s ) = 0 , which, together with (3.9), gives(3.10) E ( t (1) ) = g (0) . Therefore, the obtention of E ( t (1) ) reduces to solve the boundary value problem (3.8). In fact,(3.11) g ( x ) = − x σ for µ = 0 , − xµ − σ ( e − µxσ − e − µσ )2 µ for µ > . Moreover, one can validate that − xµ − σ ( e − µxσ − e − µσ )2 µ → − x σ as µ →
0. Hence, we conclude E ( t (1) ) = g (0) = σ for µ = 0 , µ − σ (1 − e − µ/σ )2 µ for µ > . We have derived (3.3). Next, let us turn attention to the variance
Var ( t (1) ).In view of Var ( t (1) ) = E [( t (1) ) ] − [ E ( t (1) )] , what left is to calculate E [( t (1) ) ]. From (3.9),(3.12) ( t (1) ) =( g (0)) + Z t (1) σ ∂g∂x ( φ ( s )) dW ( s ) ! + 2 g (0) Z t (1) σ ∂g∂x ( φ ( s )) dW ( s ) . Taking the expectation of (3.12), and noting that E [ R t (1) σ ∂g∂x ( φ ( s )) dW ( s )] = 0, we have (by Ito’s isom-etry)(3.13) E [( t (1) ) ] =( g (0)) + E Z t (1) σ ∂g∂x ( φ ( s )) dW ( s ) ! = ( g (0)) + E "Z t (1) σ (cid:12)(cid:12)(cid:12)(cid:12) ∂g∂x ( φ ( s )) (cid:12)(cid:12)(cid:12)(cid:12) ds , which, together with (3.10), yields(3.14) Var ( t (1) ) = E "Z t (1) σ (cid:12)(cid:12)(cid:12)(cid:12) ∂g∂x ( φ ( s )) (cid:12)(cid:12)(cid:12)(cid:12) ds . It remains to calculate the right hand side of (3.14).For any subset A in the interval [0 , A ( x ) be the characteristic function for A (i.e., 1 A ( x ) = 1for x in A , and 1 A ( x ) = 0 for otherwise). Defining the measure(3.15) π ( A ) := E "Z t (1) A ( φ ( s )) ds / E [ t (1) ] , , TATSUYA HAYASHI AND TETSUJI TOKIHIRO we rewrite (3.14) into(3.16) Var ( t (1) ) = E "Z t (1) Z σ (cid:12)(cid:12)(cid:12)(cid:12) ∂g∂x ( φ ( s )) (cid:12)(cid:12)(cid:12)(cid:12) dx ( φ ( s )) = Z σ (cid:12)(cid:12)(cid:12)(cid:12) ∂g∂x ( x ) (cid:12)(cid:12)(cid:12)(cid:12) E "Z t (1) dx ( φ ( s )) ds = E ( t (1) ) Z σ (cid:12)(cid:12)(cid:12)(cid:12) ∂g∂x (cid:12)(cid:12)(cid:12)(cid:12) π ( dx ) . Since there exists a probability density function p ( x ) satisfying(3.17) π ( A ) = Z A p ( x ) dx, π ( dx ) = p ( x ) dx. we are left with the task of finding p ( x ). In the following, we derive p ( x ) in two cases: (i) µ = 0, (ii) µ > µ = 0. Substituting g = 1 − x into (3.6), we calculate as(3.18) − g (1) − g (0) = E ( g ( φ ( t (1) − )) − g ( φ (0)))= 0 + E Z t (1) − dW ( s ) !| {z } =0 + E Z t (1) − dL ( s ) ! = − E ( L ( t (1) )) . Substituting g ( x ) = e λx into (3.6), noting that ∂g∂x = λe λx and h µ ∂∂x + σ ∂ ∂x i g = σ λ e λx , we deduce(3.19) e λ − g (1) − g (0) = E ( g ( φ ( t (1) − )) − g ( φ (0)))= λ σ E Z t (1) e λφ ( s ) ds ! + σ E Z t (1) λe λφ ( s ) dW ( s ) !| {z } =0 + E Z t (1) λe λφ ( s ) dL ( s ) ! . Since L ( t ) increases only when φ ( t ) = 0,(3.20) E Z t (1) λe λφ ( s ) dL ( s ) ! = E Z {
0. Via a similar argument to (i), instead of (3.23), one can obtain that: Z e λx ( µλ + σ λ ) p ( x ) dx = 1 E ( t (1) ) ( e λ − − λ ) + µλ, with E ( t (1) ) = µ − σ (1 − e − µ/σ )2 µ from (3.3). Then, one can verify that(3.25) p ( x ) = θ ( e θ − e θx )1 + θe θ − e θ for µ > , θ = 2 µ/σ . Hence, we have obtained the density function p ( x ) for µ ≥
0. Substituting (3.25) into (3.16), andtogether with (3.11), we obatin
Var ( t (1) ) = E [ t (1) ] Z σ (cid:12)(cid:12)(cid:12)(cid:12) ∂g∂x (cid:12)(cid:12)(cid:12)(cid:12) p ( x ) dx = σ for µ = 0 , − e − θ + 4 e − θ + 4 θe − θ + 2 θµ θ for µ > . Hence, we have proved (3.4). (cid:3)
Remark 3.5.
For any x ∈ [0 , g ( x ) of (3.11) represents the expectation of the beating interval of theoscillator with initial phase φ (0) = x . Remark 3.6.
Noting that φ is a renewal process for t ≥ t ( k ) ( k = 0 , , , · · · ), according to the renewaltheory (cf. [3, Chapter 9 (1.22) (2.25)]), p ( x ) of (3.17) is indeed the probability density of the distribution φ ( t ) in [0 ,
1] as t → ∞ . Let ˜ p ( x, t ) denote the probability density of the distribution of φ at time t . ˜ p ( x, t )satisfies the Fokker-Planck equation, or called the forward equation: ∂ ˜ p∂t − (cid:20) σ ∂ ∂x − µ ∂∂x (cid:21) ˜ p = 0 ( t, x ) ∈ (0 , ∞ ) × (0 , , (3.26a) ˜ p (1 , t ) = 0 , (3.26b) (cid:20) σ ∂ ˜ p∂x ( x, t ) − µ ˜ p ( x, t ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) x =1 x =0 = 0 , (3.26c) ˜ p ( x,
0) = δ ( x ) , (3.26d)where δ ( x ) denotes the Dirac Delta function and (3.26d) follows from the initial state of φ , i.e., φ (0) = 0.Since φ ( t ) jumps to 0 immediately when approaching 1, the density of φ ( t ) at x = 1 is zero and the fluxof the density at x = 0 , R ˜ p ( x, t ) dx = R ˜ p ( x, dx = 1 forall t >
0. The obtention of (3.26) follows from the classical argument (cf. [25, § t → ∞ , one can validate that ˜ p ( x, t ) converges to the stationary state, i.e., the solution of − (cid:20) − µ ∂∂x + σ ∂ ∂x (cid:21) p = 0 x ∈ (0 , , (3.27a) p (1) = 0 , (cid:20) σ ∂p∂x − µp (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) x =1 x =0 = 0 , (3.27b) Z p ( x ) dx = 1 . (3.27c)One can validate that p ( x ) given by (3.24) and (3.25) indeed satisfies (3.27) for µ = 0 and µ > , TATSUYA HAYASHI AND TETSUJI TOKIHIRO The phase models for two coupled cells
As explained in Section 2, the Kuramoto model is an applicable tool to investigate the synchronizationbeating of two-coupled cardiomyocytes. The conventional Kuramoto model with noise effect for tow-coupled oscillators { ¯ φ i } i =1 , is presented as follows: d ¯ φ i ( t ) = µ i dt + A i,j f ( ¯ φ j − ¯ φ i ) dt + σ i dW i ( t ) , (4.1a) ¯ φ i (0) = 0 , (4.1b)where i, j = 1 , i = j , f ( x ) = sin(2 πx ), ( µ i , σ i ) denotes the intrinsic frequency and noise strength for cell(oscillator) i , A i,j the coefficient describing the strength of reaction between cell i and cell j ( A i,j ≥ { W i ( t ) } i =1 , the two independent standard Brownian motions.However, in general case, the above model may be inadequate to capture the essential propertiesof cardiomyocytes’ synchronization. First, the irreversibility of beating should be taken into account.Second, the cardiomyocyte can be induced to beat by the neighboring cells’ action potential. In addition,after beating the cardiomyocyte enters into a refractory , during which the cell cannot be induced to beat.The length of refractory depends on the membrane potential, or more precisely, the concentrations ofCa , K + , Na + ions interior and exterior of the membrane.To incorporate the irreversibility of beating, induced beating and refractory, we modify the conventionalKuramoto model (2.11) as follows. Let { φ i } i =1 , be the phase of two cardiomyocytes, satisfying dφ i ( t ) = µ i dt + A i,j f ( φ j − φ i ) dt + σ i dW i ( t ) + dL i ( t ) , (4.2a) φ i (0) = 0 , (4.2b)where the process L i ( t ) imposes the reflective boundary for φ i . Let t ( k ) i be the k -th passage time that cell i beats. Then we call [ t ( k ) i , t ( k +1) i ) the k -th beating interval (or oscillation cycle) of cell i . L i is defined by:(L1) L i is continuous and nondecreasing during each beating interval of φ i ;(L2) L i increases only when φ i = 0.Since the refractory period associates with the membrane potential ( u of FN model), which correspondsto the phase φ i , for simplicity, we set a refractory threshold B i (0 ≤ B i < i is out of refractory, cell i beats promptly when the neighbor (cell j ) beats spontaniously,in other words, if φ i ( t − ) > B i and φ j ( t − ) = 1, then φ j ( t ) = φ i ( t ) = 0 (both two phases jump to0 after beating to start a new oscillation cycle);(REF) If cell i is in refractory and the neighbor cell j beats spontaniously, then cell i will not be inducedto beat, namely, if φ i ( t − ) ≤ B i and φ j ( t − ) = 1, then φ j ( t ) = 0 and φ i ( t ) = φ i ( t − ) ( φ j jumps to0 but φ i keeps going). Remark 4.1.
One weak point of the conventional model (4.1) is that the synchronization has been treated“ambiguous” or “approximately”, because the possibility of ¯ φ i ( t ) = ¯ φ j ( t ) = 1 is zero, and one can onlyexpect that both two cells beat with tiny time-delay, namely, ¯ φ i ( t i ) = 1, ¯ φ j ( t j ) = 1 with | t i − t j | ≈
0. Toguarantee this “approximated” synchronization, one should take sufficiently small noise strength σ i andlarge enough reaction coefficient A i,j (see Section 4.2).Thanks to the induced beating (IND), we have a rigorous mathematical definition of the synchroniza-tion. Let ≈ ( k ) denote the time of k -th synchronized beating, i.e.,(4.3) ≈ ( k ) := inf { t > ≈ ( k − : φ i ( t − ) = 1 , φ j ( t ) > B j , i = 1 or 2 , j = i } ( ≈ (0) = 0) . In view of φ ( ≈ ( k ) ) = φ ( ≈ ( k ) ) = 0, Φ( t ) := ( φ ( t ) , φ ( t )) is a renewal stochastic process for t > ≈ ( k ) ( k = 0 , , , · · · ). Therefore, the beating intervals { ≈ ( k +1) − ≈ ( k ) } k ≥ are independent and identicallydistributed (i.i.d.). To obtain the expected value and variance of synchronized beating interval, we onlyneed to investigate ≈ (1) . STOCHASTIC PHASE MODEL FOR CARDIAC MUSCLE CELLS 13
In Figure 4.1 (a)(b) and (c)(d), we plot two examples of ( ¯ φ , ¯ φ ) and ( φ , φ ) respectively. For smallnoise strength and large reaction coefficients, for example, σ = σ = 0 . A , = A , = 6, the conventionalmodel (4.1) (Figure 4.1 (b)) and the proposed model (4.2) (L1)(L2)(IND)(REF) (Figure 4.1 (d)) havesimilar solution behavor. In view of Figure 4.1 (b), the non-positive phase ( ¯ φ i ≤
0) is ignorable, and thetime-delay between two cells’ beating is very tiny. Therefore, the roles of the reflective boundary andinduced beating of our model are negligible. However, when the noise strength is not so small and thereaction coefficients is not large enough, the conventional model ( ¯ φ , ¯ φ ) may have no synchronization(see Figure 4.1 (a)). wherea Figure 4.1 (c) shows the synchronization owing to the induced beating, andthe significant role of the reflective boundary. (a) (b)(c) (d) Figure 4.1.
The trajectories of ( ¯ φ , ¯ φ ) and ( φ , φ ) with ( µ , µ ) = (1 , A , = A , = A , σ = σ = σ , B = B = B = 0 .
3. (a)(c) ( σ, A ) = (1 , σ, A ) = (0 . , ≈ (1) . First,for the proposed model, applying the Ito’s calculus and the renewal theory, we derive the PDEs associatedwith E ( ≈ (1) ) and Var ( ≈ (1) ). However, the closed-form of the PDEs’ solutions are non-trivial to derive.Next, we consider the case that the role of the reflective boundary and induced beating is negligible, wherethe conventional and proposed models have little difference, and we obtain the relationship between theparameters and the CV of beating intervals.4.1. The expectation and variance of the synchronized beating interval.
Besides of the syn-chronized beating (induced beating (IND)), let us pay attention to the single-beating (REF), or calledindependent-beating, where only one cell is beating and the other is in refractory. Assume that beforethe first synchronized beating, cell i beats independently for J i times. Let t ( m ) i ( m ≤ J i ) be the passagetime of m -th single-beating of cell i , that is,(4.4) t ( m ) i := inf { t > t ( m − i : φ i ( t − ) = 1 , φ j ( t ) < B j , t < ≈ (1) } , t (0) i = 0 ( j = i ) . For t < ≈ (1) , φ i ( t ) is a stochastic process with jump at { t ( m ) i } J i m =1 , where φ i ( t ( m ) i − ) = 1 and φ i ( t ( m ) i ) = 0.Setting Φ( t ) := ( φ ( t ) , φ ( t )), and applying the Itˆo’s formula for the stochastic process with jumps (cf. , TATSUYA HAYASHI AND TETSUJI TOKIHIRO [12]), we have: for any g ( x , x ) (0 ≤ x , x ≤
1) with continuous differential ∂g∂x i and ∂ g∂x i ∂x j ( i, j = 1 , g (Φ( ≈ (1) − )) = g (Φ(0)) + Z ≈ (1) X i =1 (cid:20) ( µ i + A i,j f ( x j − x i )) ∂∂x i + σ i ∂ ∂x i (cid:21) g (Φ( s )) ds + Z ≈ (1) X i =1 σ i ∂g∂x i (Φ( s )) dW i ( s ) + X i =1 Z ≈ (1) ∂g∂x i (Φ( s )) dL i ( s ) − J X m =1 h g (1 , φ ( t ( m )1 − )) − g (0 , φ ( t ( m )1 )) i − J X m =1 h g ( φ ( t ( m )2 − ) , − g ( φ ( t ( m )2 ) , i , where φ i ( t ( m ) j ) = φ i ( t ( m ) j − ) ≤ B i , that is, cell j is in refractory while cell i is beating ( φ i ( t ( m ) i − ) = 1).Since L i ( t ) is non-decreasing in the time intervals { ( t ( m ) i , t ( m +1) i ) } J i − m =0 and ( t ( J i ) i , ≈ (1) ), and L i ( t ) in-creases only when φ i ( t ) = 0, we see that(4.6) Z ≈ (1) ∂g∂x i (Φ( s )) dL i ( s ) = Z {
1) = 0 for B < x ≤ , (4.7d) g ( x ,
0) = g ( x ,
1) for 0 ≤ x ≤ B , (4.7e) g (0 , x ) = g (1 , x ) for 0 ≤ x ≤ B . (4.7f)It follows from (4.7b) and (4.6) that Z ≈ (1) ∂g∂x i (Φ( s )) dL i ( s ) = 0 . Since the induced beating happens at t = ≈ (1) , we have φ ( ≈ (1) − ) = 1 , φ ( ≈ (1) − ) > B or φ ( ≈ (1) − ) =1 , φ ( ≈ (1) − ) > B , which, together with (4.7c) and (4.7d), implies g (Φ( ≈ (1) − )) = 0 . Furthermore, in view of φ i ( t ( m ) j ) = φ i ( t ( m ) j − ) ≤ B i , (4.7e) and (4.7f) guarantee that J X m =1 h g (1 , φ ( t ( m )1 − )) − g (0 , φ ( t ( m )1 )) i = 0 , J X m =1 h g ( φ ( t ( m )2 − ) , − g ( φ ( t ( m )2 ) , i = 0 . With Φ(0) = (0 ,
0) and (4.7a), we rewrite (4.5) into:(4.8) 0 = g (0 , − ≈ (1) + Z ≈ (1) X i =1 σ i ∂g∂x i (Φ( s )) dW i ( s )Taking the expectation, we have(4.9) E [ ≈ (1) ] = g (0 ,
0) + E "Z ≈ (1) X i =1 σ i ∂g∂x i (Φ( s )) dW i ( s ) =0 = g (0 , . STOCHASTIC PHASE MODEL FOR CARDIAC MUSCLE CELLS 15
Hence, the obtention of E ( ≈ (1) ) reduces to solve the PDE (4.7) Next, let us turn attention to the variance Var [ ≈ (1) ]. From (4.8), we have( ≈ (1) ) = ( g (0 , + "Z ≈ (1) X i =1 σ i ∂g∂x i (Φ( s )) dW i ( s ) + 2 g (0 , Z ≈ (1) X i =1 σ i ∂g∂x i (Φ( s )) dW i ( s ) . Taking the expectation of the above equation yields: E [( ≈ (1) ) ] = ( g (0 , + E "Z ≈ (1) X i =1 σ i ∂g∂x i (Φ( s )) dW i ( s ) + 2 g (0 , E "Z ≈ (1) X i =1 σ i ∂g∂x i (Φ( s )) dW i ( s ) =0 = ( g (0 , + E "Z ≈ (1) X i =1 σ i (cid:12)(cid:12)(cid:12)(cid:12) ∂g∂x i (Φ( s )) (cid:12)(cid:12)(cid:12)(cid:12) ds (by Ito’s isometry) , which, together with Var [ ≈ (1) ] = E [( ≈ (1) ) ] − ( E [( ≈ (1) )]) and (4.9), implies(4.10) Var [ ≈ (1) ] = E "Z ≈ (1) X i =1 σ i (cid:12)(cid:12)(cid:12)(cid:12) ∂g∂x i (Φ( s )) (cid:12)(cid:12)(cid:12)(cid:12) ds . Noting that Φ( t ) = ( φ ( t ) , φ ( t )) returns to (0 ,
0) after every synchronization, Φ( t ) is a renewal process.According to the renewal theory [3], the right-hand side of (4.10) is evaluated as: Var ( ≈ (1) ) = E ( ≈ (1) ) Z Z
10 2 X i =1 σ i (cid:12)(cid:12)(cid:12)(cid:12) ∂g∂x i (cid:12)(cid:12)(cid:12)(cid:12) p dx dx , where p ( x , x ) denotes the distribution density of Φ = ( φ , φ ) ∈ [0 , as t → ∞ .To derive the equations concerned about p ( x , x ), we interpret the model (4.2) (L1)(L2)(IND)(REF)from a physical point of view. Regard ( φ , φ ) as the position of a particle in [0 , , which moves with ve-locity [ µ + A , f ( φ ( t ) − φ ( t )) , µ + A , f ( φ ( t ) − φ ( t ))] ⊤ , and is effected by noise [ σ dW ( t ) , σ dW ( t )] ⊤ .The initial position of the particle is (0 , { x = 1 , x > B }∪{ x >B , x = 1 } , then it jumps to point (0 ,
0) immediately. On the other hand, if the particle approaches theboundary { x = 1 , < x ≤ B } (resp. { < x ≤ B , x = 1 } ), it jumps to position (0 , x ) (resp. ( x , { ( x , x ) | x = 0 or x = 0 } . InFigure 4.2 (a)(b), we plot two trajectories of the particle.Therefore, p ( x , x ) represents the distribution density of the particle in [0 , as t → ∞ . Now, let˜ p ( x , x , t ) denote the distribution density of the particle at time t . Via a similar argument to [25, seeSection 3.5], one can prove that ˜ p ( x , x , t ) satisfies the following Fokker-Planck equation (or the forwardequation): for i, j = 1 , i = j , ∂ ˜ p∂t + X i =1 ∂∂x i ˜ F i = δ ( x , x ) X i =1 Z { x i =1 , B j ≤ x j < } ˜ F i ds for 0 < x , x < , t > , (4.11a) ˜ p ( x , x , t ) = 0 on { x = 1 } ∪ { x = 1 } , (4.11b) ˜ F = 0 on { x = 0 , B ≤ x < } , (4.11c) ˜ F = 0 on { x = 0 , B ≤ x < } , (4.11d) ˜ F | x =1 x =0 = 0 for { < x ≤ B } , (4.11e) ˜ F | x =1 x =0 = 0 for { < x ≤ B } , (4.11f) ˜ p ( x , x ,
0) = δ ( x , x ) , for { < x , x < } , (4.11g)where ˜ F i := ( µ i + A i,j f ( x j − x i ))˜ p − σ i ∂ ˜ p∂x i denotes the i -th component of flux, and δ ( x , x ) the DiracDelta function. (4.11g) means the initial position of the particle is (0 , , TATSUYA HAYASHI AND TETSUJI TOKIHIRO (a) (b) Figure 4.2.
Two realizations of the trajectories ( φ , φ ) with B = B = 0 .
3. Thereflective boundary is imposed on { x = 0 } ∪ { x = 0 } . (a) The particle touches theboundary (blue color) { x > B , x = 1 } , which will jumps to point (0 , { < x ≤ B , x = 1 } (cyan trajectory), and jumpsto point ( x , { x > B , x = 1 } (red trajectory), whichwill jumps to point (0 , φ , φ ) correspond to the induced bearing(IND), whereas the cyan one corresponds to the independent beating (REF).(0 ,
0) or ∪ i = ji,j =1 , { x i = 0 , < x j < B j } instantly when touching the boundary ∪ i =1 , { x i = 1 } , the densityof the particle on ∪ i =1 , { x i = 1 } is zero, namely (4.11b), which implies X i =1 Z { x i =1 ,B j ≤ x j < } − σ i ∂ ˜ p∂x i ds = X i =1 Z { x i =1 ,B j ≤ x j < } ˜ F i ds. Hence, the right-hand side of (4.11a) represents the total flux of ˜ p which touches the boundary { x =1 , B ≤ x < }∪{ B ≤ x < , x = 1 } and then jumps to point (0 ,
0) immediately. Putting together withthe boundary conditions (4.11c)–(4.11g), and in view of R [0 , ˜ p ( x , x , dx = R [0 , δ ( x , x ) dx = 1,one can validate the conservation law: ddt Z [0 , ˜ p ( x , x , t ) dx = 0 (equivalently Z [0 , ˜ p ( x , x , t ) dx = 1) for all t > . The zero-flux boundary conditions (4.11c), (4.11d) correspond to the reflective boundary, whereas (4.11e)(resp. (4.11f)) describes that the particle jumps from (1 , x ) to (0 , x ) with 0 < x ≤ B (resp. from( x , x ) = ( x ,
1) to ( x ,
0) with 0 < x ≤ B ).As t → ∞ , one can show that ˜ p ( x , x , t ) converges to the stationary state p ( x , x ), which satisfies:for i, j = 1 , i = j , X i =1 ∂∂x i F i = δ ( x , x ) X i =1 Z { x i =1 , B j ≤ x j < } F i ds for 0 < x , x < , (4.12a) p ( x , x ) = 0 on { x = 1 } ∪ { x = 1 } , (4.12b) F = 0 on { x = 0 , B ≤ x < } , (4.12c) F = 0 on { x = 0 , B ≤ x < } , (4.12d) F | x =1 x =0 = 0 for { < x ≤ B } , (4.12e) F | x =1 x =0 = 0 for { < x ≤ B } , (4.12f) STOCHASTIC PHASE MODEL FOR CARDIAC MUSCLE CELLS 17 where F i := ( µ i + A i,j f ( x j − x i )) p − σ i ∂p∂x i denotes the i -th component of flux ( F i = − σ i ∂p∂x i on { x i = 1 } by (4.12b)). As the stationary state of ˜ p , p also satisfies R [0 , p dx dx = 1.From the above argument, we conclude: Proposition 4.1.
The expectation and variance of the synchronized beating interval are given by E ( ≈ (1) ) = g (0 , , (4.13a) Var ( ≈ (1) ) = E ( ≈ (1) ) Z [0 , X i =1 σ i (cid:12)(cid:12)(cid:12)(cid:12) ∂g∂x i (cid:12)(cid:12)(cid:12)(cid:12) p dx dx , (4.13b) where g, p are the solutions to (4.7) and (4.12) , respectively. In the case of single isolated cardiomyocyte (Section 3), we have obtained g and p in closed-form.However, it is non-trivial to solve the two-dimensional PDEs (4.7) and (4.11). On the elementary casethat µ = µ = 0, A , = A , = 0 and B = B = 0, (4.7) is reduced to: X i =1 σ i ∂ g∂x i = − < x , x < , (4.14a) ∂g∂x i = 0 for x i = 0 , i = 1 , , (4.14b) g (1 , x ) = 0 for 0 < x ≤ , (4.14c) g ( x ,
1) = 0 for 0 < x ≤ . (4.14d)Apparently, the eigenvalues { λ mn } ∞ m,n =0 and eigenfunctions { g mn } ∞ m,n =0 for theoperator - P i =1 σ i ∂ ∂x i un-der the boundary conditions (4.14b)–(4.14d) are given by: λ mn = σ m + 12 ) π + σ n + 12 ) π , g mn = cos(( m + 12 ) πx ) cos(( n + 12 ) πx ) . Then, there exist constants { a mn } ∞ m,n =0 such that g = P ∞ m,n =0 a mn g mn is the solution of (4.14). Substi-tuting g = P ∞ m,n =0 a mn g mn into (4.14a), and calculating the integration R [0 , (4.14a) × g m ′ n ′ dx dx for m ′ , n ′ = 0 , , , · · · , one can derive that a mn = 4( − m + n λ mn ( m + )( n + ) π , m, n = 0 , , , · · · . Hence, we get the expected value of the synchronized beating interval E ( ≈ (1) ) = g (0 ,
0) = ∞ X m,n =0 a mn g mn (0 ,
0) = ∞ X m,n =0 − m + n λ mn ( m + )( n + ) π < ∞ . In above, we have derive E ( ≈ (1) ) for the case with zero intrinsic frequencies { µ i } , zero reaction coefficients { A i,j } , and zero refractory thresholds { B i } . However, for the general case, the closed-form of g and p are difficult to obtain, where one can compute the numerical solutions using the finite difference/elementmethod (see Figure 4.3 for a numerical example of g and p ).4.2. The synchronized beating of the conventional model.
In view of Figure 4.1 (b)(d), when thenoise strength is sufficiently small and the reaction coefficients are large enough, the role of the reflectiveboundary and induced beating is ignorable, such that there is no much difference between the proposedmodel (4.2) (L1)(L2)(IND)(REF) and the conventional model (4.1).In this section, we shall pay attention to the conventional model (4.1). Since the probability for the“exact synxhronization” ¯ φ ( t ) = ¯ φ ( t ) = 1 is zero, we can only consider the “approximated synchroniza-tion”, i.e., ¯ φ ( t ) = ¯ φ ( t ) = 1 with t ≈ t . Let the k -th beating time for oscillator i be the k -th passagetime that ¯ φ i = 1: t ( k ) i := inf { t > t ( k − i : ¯ φ i ( t ) = 1 } ( t (0) i = 0) . , TATSUYA HAYASHI AND TETSUJI TOKIHIRO (a) (b) Figure 4.3. (a) The profile of g ( x , x ) with E ( ≈ (1) ) = g (0 , ≈ .
43. (b) The profileof p ( x , x ). Here, we set ( µ , µ ) = (1 , σ = σ = 1, A , = A , = 2, B = B = 0 . Remark 4.2.
In view of f ( ¯ φ j − ¯ φ i ) = sin(2 π ( ¯ φ j − ¯ φ i )), it is equivalent that we remove the setting that¯ φ i jumps to 0 when reaching 1, and define the k -th “synchronized” beating time t ( k ) i as the first passagetime that ¯ φ i = k . For the convenience of the discussion, we temporarily remove the enforcement that¯ φ i ( t ) = 0 if ¯ φ i ( t − ) = 1 in the following argument of this section. Hence, the k -th beating time of oscillator i is redefined by: t ( k ) i := inf { t > t ( k − i : ¯ φ i ( t ) = k } . In addition, we assume the “approximated synchronization” occurs, saying t ( k )1 ≈ t ( k )2 .To ensure the “approximated synchronization”, we assume that | ¯ φ − ¯ φ | ≪
1. In fact, we show thatfor sufficiently large reaction coefficients { A i,j } and small enough noise strength { σ i } , one can guaranteethat | E ( ¯ φ − ¯ φ ) | ≤ ǫ ≪ Var ( ¯ φ − ¯ φ ) ≤ ǫ ≪ d ¯ φ = µ dt + A , sin(2 π ( ¯ φ − ¯ φ )) dt + σ dW ( t ) ,d ¯ φ = µ dt + A , sin(2 π ( ¯ φ − ¯ φ )) dt + σ dW ( t ) , we get d ( ¯ φ − ¯ φ ) = ( µ − µ ) dt + ( A , + A , ) sin(2 π ( ¯ φ − ¯ φ )) dt + σ dW ( t ) − σ dW ( t ) . For | ¯ φ − ¯ φ | ≪
1, we adopt the approximation sin(2 π ( ¯ φ ( t ) − ¯ φ ( t ))) ≈ π ( ¯ φ ( t ) − ¯ φ ( t )). Then the aboveequation becomes d [( ¯ φ − ¯ φ ( t ))] = ( µ − µ ) dt + 2 π ( A , + A , )( ¯ φ − ¯ φ ) dt + σ dW ( t ) − σ dW ( t ) , which is equivalent to d [ e π ( A , + A , ) t ( ¯ φ − ¯ φ ) = e π ( A , + A , ) t [( µ − µ ) dt + σ dW ( t ) − σ dW ( t )] . With the initial value ¯ φ (0) − ¯ φ (0) = 0, we find that(4.15) ¯ φ ( t ) − ¯ φ ( t ) = ( µ − µ ) 1 − e − π ( A , + A , ) t π ( A , + A , ) + Z t e π ( A , + A , )( s − t ) [ σ dW ( s ) − σ dW ( s )] . Taking the expectation of (4.15) yields(4.16) E [ ¯ φ ( t ) − ¯ φ ( t )] = ( µ − µ ) 1 − e − π ( A , + A , ) t π ( A , + A , ) + 0 . Thus, for sufficiently large { A i,j } such that 2 π ( A , + A , ) ≥ | µ − µ | ε (0 < ε ≪ | E [ ¯ φ − ¯ φ ] | ≤ ε isguaranteed. STOCHASTIC PHASE MODEL FOR CARDIAC MUSCLE CELLS 19
To derive the sufficient condition for
Var ( ¯ φ − ¯ φ ) ≤ ǫ ≪
1, from (4.15), (4.16), we calculate as( ¯ φ ( t ) − ¯ φ ( t )) =( E [ ¯ φ ( t ) − ¯ φ ( t )]) + (cid:18)Z t e π ( A , + A , )( s − t ) [ σ dW ( s ) − σ dW ( s )] (cid:19) + 2 E [ ¯ φ ( t ) − ¯ φ ( t )] Z t e π ( A , + A , )( s − t ) [ σ dW ( s ) − σ dW ( s )] , which implies E [( ¯ φ − ¯ φ ) ] = ( E [ ¯ φ − ¯ φ ]) + E (cid:20)Z t e π ( A , + A , )( s − t ) σ dW ( s ) (cid:21) + E (cid:20)Z t e π ( A , + A , )( s − t ) σ dW ( s ) (cid:21) + 2 E (cid:20)Z t e π ( A , + A , )( s − t ) σ dW ( s ) (cid:21)| {z } =0 E (cid:20)Z t e π ( A , + A , )( s − t ) σ dW ( s ) (cid:21)| {z } =0 + 2 E [ ¯ φ ( t ) − ¯ φ ( t )] E (cid:20)Z t e π ( A , + A , )( s − t ) [ σ dW ( s ) − σ dW ( s )] (cid:21)| {z } =0 By Ito’s isometry, we have
Var [ ¯ φ − ¯ φ ] = E [( ¯ φ − ¯ φ ) ] − ( E [ ¯ φ − ¯ φ ]) = Z t e π ( A , + A , )( s − t ) ( σ + σ ) ds = (1 − e − π ( A , + A , ) t )( σ + σ ) . Therefore, for sufficiently small noise strength such that ( σ + σ ) ≤ ε , we have Var [ ¯ φ ( t ) − ¯ φ ( t )] ≤ ǫ .From now on, we tacitly assume that { A i,j } are sufficiently large and { σ i } are small enough suchthat the “approximated” synchronization ( t ( k )1 ≈ t ( k ) j ) occurs. And we turn to investigate the CV of thebeating intervals { t ( k ) i } , where we employ the approximation approach proposed by [21, (5)–(18)].Let us briefly introduce the idea of [21]. For a very large time scale, one can approximate the stablesynchronization oscillation system by the linear system:(4.17) φ syn i ( t ) = µ syn t + ψ syn i , i = 1 , , where the phase functions { φ syn i } i =1 , are called the synchronized solutions, with the intrinsic synchronizedfrequency µ syn and initial state ψ syn i satisfying(4.18) A i,j sin(2 π ( ψ syn j − ψ syn i )) = µ i − µ syn , i, j = 1 , , i = j. For (4.17), we have the synchronized beating interval τ = 1 /µ syn , and(4.19) 1 = φ syn i ( t + τ ) − φ syn i ( t ) . Here, we take τ as the mean value of the beating intervals { t ( k ) i } k =1 , , ··· . Since one oscillation cycle of ¯ φ i corresponds to the increasement of ¯ φ i by 1, according to the discussion of [21], the variance of beatingintervals { t ( k ) i } k =1 , , ··· is proportional to the variance of ¯ φ i ( t + τ ) − ¯ φ i ( t ) − t → ∞ . Therefore, theCV of { t ( k ) i } can be approximated by(4.20) CV i := q lim t →∞ E [( ¯ φ i ( t + τ ) − ¯ φ i ( t ) − ] . Setting the notation ξ i ( t ) := ¯ φ i ( t ) − φ syn i ( t ), from (4.19) and (4.20), we see that(4.21) CV i = lim t →∞ E [( ξ i ( t + τ ) − ξ i ( t )) ] , TATSUYA HAYASHI AND TETSUJI TOKIHIRO Now the problem reduces to calculate E [( ξ i ( t + τ ) − ξ i ( t )) ]. To this end, we first derive the equationsfor { ξ i } i =1 , : i, j = 1 , i = j , dξ i ( t ) = ( µ i − µ syn ) dt + A i,j sin(2 π ( ¯ φ j ( t ) − ¯ φ i ( t ))) dt + σ i dW i ( t ) ,ξ i (0) = ξ i := ¯ φ i (0) − φ syn i (0) = − ψ syn i , We assume that the difference between the synchronized solution φ syn i and phase ¯ φ i is small, i.e., | ξ i | = | φ syn i − ¯ φ i | ≪ . In view ofsin(2 π ( ¯ φ j ( t ) − ¯ φ i ( t ))) = sin(2 π ( ξ j ( t ) − ξ i ( t ) + φ syn j ( t ) − φ syn i ( t ))) = sin(2 π ( ξ j ( t ) − ξ i ( t ) + ψ syn j − ψ syn i ))= sin(2 π ( ψ syn j − ψ syn i )) + cos(2 π ( ψ syn j − ψ syn i ))( ξ j ( t ) − ξ i ( t )) + O (( ξ j ( t ) − ξ i ( t )) ) , and neglecting the smaller quadratic term O ( | ξ j ( t ) | ) and O ( | ξ i ( t ) | ), together with (4.18), we obtain: dξ i ( t ) = b ij ( ξ j ( t ) − ξ i ( t )) dt + σ i dW i ( t ) , (4.22a) ξ i (0) = ξ i , (4.22b)where b ij := A i,j cos(2 π ( ψ syn j − ψ syn i )). In the following, we assume A i,j > | ψ syn j − ψ syn i | ≪ π ( ψ syn j − ψ syn i )) ≈ b ij ≈ A i,j > t →∞ E [( ξ i ( t + τ ) − ξ i ( t )) ]. Second, our result shows a explicit relationship between theparameters b ij and the CV, which is of practical use to determine the suitable parameters { A i,j } (seeRemark 4.3). Proposition 4.2.
We approximate the CV of the synchronized beating intervals { t ( k ) i } k =1 , , ··· by CV i =lim t →∞ E [( ξ i ( t + τ ) − ξ i ( t )) ] , where { ξ i } i =1 , is the solution of (4.22) . For A , , A , > , and cos(2 π ( ψ syn j − ψ syn i )) > , we have: i, j = 1 , , i = j , (4.23) CV i := lim t →∞ E [( ξ i ( t + τ ) − ξ i ( t )) ]= b − n ( b ij σ j + b ji σ i ) τ + b − (cid:2) b ij ( b ji σ i − b ij σ j ) + b ij ( σ j + σ i ) (cid:3) (1 − e − τb ) o . Remark 4.3.
In Section 3, we determine the intrinsic frequency and noise strength ( µ i , σ i ) for single-isolated cell i ( i = 1 ,
2) by formulas (3.3)–(3.5), together with the mean value and variance/CV of thebeating intervals obtained from the bio-experiments [20]. Coupling two cells (cell 1 and cell 2), we intendto find suitable coefficients A , and A , for the reaction terms. Assuming that the difference between thesynchronized solution { φ syn i } i =1 , is tiny ( | φ syn1 − φ syn2 | = | ψ syn1 − ψ syn2 | ≈ b ij = A i,j cos(2 π ( ψ syn i − ψ syn j )) ≈ A i,j , we see that(4.24) CV i =( A i,j + A j,i ) − ( A i,j σ j + A j,i σ i ) τ + ( A i,j + A j,i ) − (cid:2) A i,j ( A j,i σ i − A i,j σ j ) + A i,j ( σ j + σ i ) (cid:3) (1 − e − τ ( A i,j + A j,i ) ) . Meanwhile, the expetation and CV of the synchronized beating intervals, denoted by T and CV , can beobtained from the bio-experiments [20]. Substituting τ = T and CV = CV = CV into (4.24), one cansolve (4.24) ( i, j = 1 , i = j ) numerically to get the coefficients A , and A , . Proof of Proposition 4.2.
Setting the notations ξ = (cid:20) ξ ξ (cid:21) , B = (cid:20) b − b − b b (cid:21) , STOCHASTIC PHASE MODEL FOR CARDIAC MUSCLE CELLS 21 ξ = (cid:20) ξ ξ (cid:21) , W = (cid:20) W W (cid:21) , σ = (cid:20) σ σ (cid:21) , we write (4.22) as follows: d ξ = − Bξ dt + σ d W ( t ) , ξ (0) = ξ . Multiplying the above equation with e − B t , we have d ( e B t ξ ) = e B t [ ν dt + σ d W ] , which implies(4.25) ξ ( t ) = e − B t ξ + Z t e − B ( t − s ) σ d W ( s ) . Since the expectation of Itˆo’s integral is zero, E [ ξ ( τ )] = e − B t ξ . One can validate that B has two sets of eigenvalue and eigenvector: λ = 0 u = [1 , ⊤ , λ = b := b + b u = [ b , − b ] ⊤ . And we have e − t B u = u , e − t B u = e − tb u With the help of ( λ i , u i ) i =1 , , we make the decompositions ξ = b − ( b ξ + b ξ ) u + b − ( ξ − ξ ) u , σ d W = b − ( b σ dW + b σ dW ) u + b − ( σ dW − σ dW ) u , substituting which into (4.25), we observe that(4.26) ξ ( t ) = b − ( b ξ + b ξ ) + b − ( ξ − ξ ) e − tb b + Z t b − ( b σ dW ( s ) + b σ dW ( s )) + Z t b − e − ( t − s ) b b ( σ dW ( s ) − σ dW ( s ))= b − C ( t ) + b − σ Z t ( b + e − ( t − s ) b b ) dW ( s ) | {z } =: W ( t ) + b − σ Z t ( b − e − ( t − s ) b b ) dW ( s ) | {z } =: W ( t ) , where C ( t ) := ( b ξ + b ξ ) + b − ( ξ − ξ ) e − tb b . Then, we see that ξ ( t ) = b − [ C ( t ) + W ( t ) + W ( t )] + 2 b − [ C ( t ) W ( t ) + C ( t ) W ( t ) + W ( t ) W ( t )]The expectation of Itˆo’s integral W i is zero, that is, E [ W ( t )] = 0 , E [ W ( t )] = 0 , together with the independency between W ( t ) and W ( t ) (because W and W are independent), whichgives E [ W ( t ) W ( t )] = E [ W ( t )] E [ W ( t )] = 0 . Moreover, by Ito’s isometry, E [ W ( t )] = Z t ( b + e − ( t − s ) b b ) σ ds, E [ W ( t )] = Z t ( b − e − ( t − s ) b b ) σ ds. Hence, we conclude(4.27) E [ ξ ( t )] = b − C ( t ) + b − Z t ( b + e − ( t − s ) b b ) σ ds + b − Z t ( b − e − ( t − s ) b b ) σ ds = b − ( b ξ + b ξ + ( ξ − ξ ) e − tb b ) + b − (cid:20) ( b σ + b σ ) t + 2 b ( b σ − b σ ) 1 − e − tb b + b ( σ + σ ) 1 − e − tb b (cid:21) . , TATSUYA HAYASHI AND TETSUJI TOKIHIRO In view of E [( ξ ( t + τ ) − ξ ( t )) ] = E [ ξ ( t + τ )+ ξ ( t ) − ξ ( t + τ ) ξ ( t )], it remains to calculate E [ ξ ( t + τ ) ξ ( t )]. E [ ξ ( t + τ ) ξ ( t )] = b − E [ C ( t ) C ( t + τ ) + W ( t ) W ( t + τ ) + W ( t ) W ( t + τ )]+ b − E [ C ( t ) W ( t + τ ) + C ( t ) W ( t + τ ) + W ( t ) W ( t + τ ) + W ( t ) W ( t + τ )] | {z } =0 . Let us pay attention to E [ W ( t ) W ( t + τ )]. We divide W ( t + τ ) into W ( t + τ ) = σ Z t ( b + e − ( t + τ − s ) b b ) dW ( s ) + σ Z t + τt ( b + e − ( t + τ − s ) b b ) dW ( s ) . The independency between R t + τt ( b + e − ( t + τ − s ) b b ) dW ( s ) and R t ( b + e − ( t − s ) b b ) dW ( s ) yields E [ W ( t ) W ( t + τ )] = σ (cid:20)(cid:18)Z t ( b + e − ( t − s ) b b ) dW ( s ) (cid:19) (cid:18)Z t ( b + e − ( t + τ − s ) b b ) dW ( s ) (cid:19)(cid:21) = σ Z t ( b + e − ( t − s ) b b )( b + e − ( t + τ − s ) b b ) ds (by Ito’s isometry) . Treating E [ W ( t ) W ( t + τ )] in a similar way, we get(4.28) E [ ξ ( t + τ ) ξ ( t )] = b − C ( t ) C ( t + τ )+ b − σ (cid:20) b t + b b (1 + e − τb ) 1 − e − tb b + b e − τb (1 − e − tb )2 b (cid:21) + b − σ (cid:20) b t − b (1 + e − τb ) 1 − e − tb b + b e − τb (1 − e − tb )2 b (cid:21) . Following from (4.27), (4.28), we find that E [( ξ ( t + τ ) − ξ ( t )) ] = b − h b ( ξ − ξ ) e − tb (1 − e − τb ) + ( b σ + b σ ) τ i + b − (cid:20) b ( b σ − b σ ) b (1 − e − τb ) (cid:21) + b − (cid:20) b ( σ + σ )2 b (2 − e − τb + e − t + τ ) b − e − tb ) (cid:21) Passing to the limit t → ∞ and in view of b = b + b >
0, we have(4.29) ( CV ) = lim t →∞ E [( ξ ( t + τ ) − ξ ( t )) ]= b − n ( b σ + b σ ) τ + b − (cid:2) b ( b σ − b σ ) + b ( σ + σ ) (cid:3) (1 − e − τb ) o . Analogously to above argument, we can calculate CV . (cid:3) The phase model for the N -cells network Let us extend the phase models of two-coupled cells to N -cells network. Figure 5.1 (a) shows twoexamples of cell-network constructed via the on-chip cellomics technology [17, 20]. Numbering the cellsby { , , · · · , N } , we denote by N i the neighbors of cell i (see Figure 5.1 (a)). In this section, we firstintroduce the phase model for N -cells network incorporating the irreversibility of beating (reflectiveboundary), induced beating and refractory.For the case with sufficiently large reaction coefficients and small enough noise strength, the proposedmodel has similar behavior to the conventional model, and the synchronization is very stable, because theeffects of reflective boundary, induced firing and refractory is ignorable (see Figure 5.2 (b)(d)). Since themassive bio-experiments (cf. [17]) reveal that the CV of the synchronized beating intervals reduces as thenetwork size increases (in other word, the synchronization is more stable if we add more cardiamyocytesto the network), we shall investigate the network-size-dependent CV of the synchronized beating intervalsby the conventional model with a similar analysis to Section 4.2. STOCHASTIC PHASE MODEL FOR CARDIAC MUSCLE CELLS 23 (a) (b) (c)
Figure 5.1. (a) Tow examples of cell network. 4-cells network (I): N i = { } ( i = 1 , , N = { , , } . 4-cells network (II): N = { } , N = { , } , N = { , } , N = { } .(c) The size-dependent beating fluctuation. The CV of the synchronized beating intervalsdecreases as the cell number increases. (d) The log-log scale of (c), where the black straightline represents ∝ N − / . N is the number of cells in the network.5.1. The phase model of N -cells network. Let ( φ i , µ i , σ i ) denote the phase, intrinsic frequency andnoise strength of cell i ( i = 1 , . . . , N ), and A i,j the coefficient of the reaction term between cell i and j . For simplicity, we consider the network that all the cells are connected with each other, that is N i = { , , · · · , N } − { i } . Then, the equations of { φ i } Ni =1 are stated as follows: for i = 1 , , · · · , N , dφ i ( t ) = µ i dt + X j ∈N i A i,j sin(2 π ( φ j − φ i )) dt + σ i dW i ( t ) + dL i ( t ) , (5.1a) φ i (0) = 0 , (5.1b)where W i ( t ) denotes the normal Brownian motion ( { W i } Ni =1 are independent), and L i ( t ) the processimplementing the reflective boundary (see (L1)(L2) of Section 4). When φ i ( t − ) = 1, we say cell i beatsspontaneously. At the same time, the neighboring cell j ( j ∈ N i ) is induced to beat if φ j ( t − ) > B j (cell j is out of refractory), where B j ∈ [0 ,
1] denotes the refractory threshold of cell j . In this case, cell i andcell j have a synchronized beating. And after beating, both two phases jump to zero, that is, φ i ( t ) = 0and φ j ( t ) = 0. On the other hand, if φ j ( t − ) ≤ B j , we say cell j is in refractory and cannot be induced tobeat, and we have φ j ( t ) = φ j ( t − ).As with Section 4, we also pay attention to the conventional Kuramoto model. Let ¯ φ i denote the phaseof cell i ( i = 1 , , · · · , N ) without the reflective boundary and induced beating, which satisfies d ¯ φ i ( t ) = µ i dt + X j ∈N i A i,j sin(2 π ( ¯ φ j − ¯ φ i )) dt + σ i dW i ( t ) , (5.2a) ¯ φ i (0) = 0 . (5.2b)In Figure 5.2 (a)(b) and (c)(d), we plot two trajectories of { φ i ( t ) } i =1 and { ¯ φ i ( t ) } i =1 respectively fordifferent parameters. For σ i = 1, A i,j = 1, the proposed model (5.1) has the synchornization due tothe induced beating (see Figure 5.2(a)), and the noise effect at φ i = 0 has been inhibited by reflectiveboundary (the irreversibility of beating). However, there is no synchronization for the conventional model(5.2) (see Figure 5.2(c)).As discussed in Section 4.2, we have to choose large enough reaction coefficients and sufficiently smallnoise strength to expect the “approximated” synchronization occurs for the conventional model. For σ i = 0 . A i,j = 8, our model (5.1) gives a very stable synchronization (see Figure 5.2(b)), where the roleof the reflective boundary and induced beating can be negligible, such that the solution behaviors of (5.1)and the conventional model (5.2) (see Figure 5.2(d)) are quite similar.In bio-experiments, the fluctuation of the synchronized beating intervals reduces as the network sizeincreases. Since both two models (5.1) and (5.2) are quite similar when the stable synchronization happens , TATSUYA HAYASHI AND TETSUJI TOKIHIRO (a) (b)(c) (d) Figure 5.2. (a)(b): The realizations of { φ i } i =1 . (c)(d): The realizations of { ¯ φ i } i =1 . Theintrinsic frequency ( µ , µ , µ , µ ) = (2 , . , , . B i = 0 . i = 1 , , ,
4) (see the dot lines in (a)(b)). (a)(c): σ i = 1, A i,j = 1; (b)(d): σ i = 0 . A i,j = 8 ( i = 1 , , , j ∈ N i ).(Figure 5.2(b)(d)), from now on, we shall pay attention to the conventional model (5.2), and derive theCV of beating interval via a similar methodology to Section 4.2.We assume that all the cells beat almost simultaneously and ignore the tiny difference between thebeating time of each ¯ φ i . Then, analogously to (4.17), taking the expected beating interval τ as thesynchronized beating interval, we introduce the synchronized solution { φ syn i } Ni =1 :(5.3) φ syn i ( t ) = µ syn t + ψ syn i , i = 1 , , · · · , N, where µ syn := 1 /τ represents the intrinsic frequency of synchronization, and ψ syn i satisfies X j ∈N i A i,j sin(2 π ( ψ syn j − ψ syn i )) = µ i − µ syn , i = 1 , , · · · , N. For the conventional model, because the reaction term is sin(2 π ( ¯ φ j − ¯ φ i )) it makes no difference that weremove the setting “the phase jumps to 0 when approaching to 1” and define the k -th beating time of cell i as the passage time that ¯ φ i reaches k (see Remark 4.2).We consider the case that σ i = σ , A i,j = A j,i for all i, j ∈ { , , · · · , N } . Via a similar approach to(4.20)–(4.22), we shall calculate(5.4) ( CV ) := 1 N N X i =1 ( CV i ) = 1 N N X i =1 lim t →∞ E [( ξ i ( t + τ ) − ξ i ( t )) ] = 1 N lim t →∞ E " N X i =1 | ξ i ( t + τ ) − ξ i ( t ) | , where ξ i := ¯ φ i − φ syn i satisfies dξ i ( t ) = X j ∈N i b ij ( ξ j ( t ) − ξ i ( t )) dt + σ i dW i ( t ) , (5.5a) ξ i (0) = ξ i , (5.5b)with b ij := A i,j cos(2 π ( ψ syn j − ψ syn i )) and ξ i = − ψ syn i . Here, we assume that A i,j > | φ syn j − φ syn i | = | ψ syn j − ψ syn i | ≪
1, such that b ij ≈ A i,j > STOCHASTIC PHASE MODEL FOR CARDIAC MUSCLE CELLS 25
Proposition 5.1.
For identical noise strength σ i = σ ( i = 1 , , · · · , N ) and symmetry reaction coefficients A i,j = A j,i , the fluctuation of the synchronized beating interval is given by: (5.6) CV √ τ = 1 √ N σ vuut N X i =2 − e − τλ i τ λ i , where CV is defined by (5.4) with ξ i satisfies (5.5) ( i = 1 , , · · · , N ). Remark 5.1.
It is known that lim N →∞ N N X i =2 − e − τλ i τ λ i = λ ∞ , where λ ∞ is some constant. Therefore, the fluctuation CV / √ τ decreases with order O ( N − ) when N is not so large, and converges to the constant σλ ∞ as N → ∞ , which has been confirmed by numericalsimulation (see Figure 5.1 (b)(c)). Proposition 5.1 is similar to the result of [21]. However, we emphasizethat we establish a new analysis with more rigorous and precious mathematical argument using thestochastic calculus. Proof of Proposition 5.1.
Setting the notations ξ = [ ξ i ] , ξ = [ ξ i ] , W = [ W i ] , B = [ − b ij ] , σ = diag( σ , · · · , σ N )(here b ii = − P j ∈N i b ij and b ij = 0 if j / ∈ N i ), we see that(5.7) d ξ = − Bξ dt + σ d W ( t ) , ξ (0) = ξ , which implies ξ ( t ) = e − t B ξ + Z t e − ( t − s ) B σ d W ( s ) . Setting | ξ | = P Ni =1 ξ i , from ξ ( t + τ ) − ξ ( t ) = [ e − ( t + τ ) B − e − t B ] ξ + Z t [ e − ( t + τ − s ) B − e − ( t − s ) B ] σ d W ( s ) + Z t + τt e − ( t + τ − s ) B σ d W , we obtain(5.8) E [ | ξ i ( t + τ ) − ξ i ( t ) | ] = | e − t B ( e − τ B − I ) ξ | + Z t | e − ( t − s ) B ( e − τ B − I ) σ | ds + Z t + τt | e − ( t + τ − s ) B σ | ds, where we have used the fact that the expectation of Itˆo’s integral is zero and { W i } are independentBrownian motion. For σ i = σ and A i,j = A j,i > i, j = 1 , , · · · , N ), σ = σ I and B = [ − b ij ] issymmetry, as well as e − t B and e − t B − I . Hence, we calculate as(5.9) | e − ( t − s ) B ( e − τ B − I ) σ | = N X i,j =1 | [ e − ( t − s ) B ( e − τ B − I ) σ ] ij | = tr( e − t − s ) B ( e − τ B − I ) σ )= N X i =1 σ ( e − τλ i − e − t − s ) λ i , (5.10) | e − ( t + τ − s ) B σ | = N X i,j =1 | [ e − ( t + τ − s ) B σ ] ij | = tr( e − t + τ − s ) B σ ) = N X i =1 σ e − t + τ − s ) λ i , where [ C ] ij denotes the ( i, j ) component of matrix C , and { λ i } Ni =1 the eigenvalues of B . , TATSUYA HAYASHI AND TETSUJI TOKIHIRO In view of b ii = − P j = i b ij , b ij > j ∈ N i and b ij = 0 for j / ∈ N i , one can validate that theeigenvalues of B = [ − b ij ] satisfies:(5.11) λ N ≥ λ N − ≥ · · · ≥ λ > λ = 0 . Substituting (5.9), (5.10) into (5.8), we obtain(5.12) E [ | ξ i ( t + τ ) − ξ i ( t ) | ] = | e − t B ( e − τ B − I ) ξ | + N X i =1 σ ( e − τλ i − Z t e − t − s ) λ i ds + N X i =1 σ e − t + τ − s ) λ i ds = | e − t B ( e − τ B − I ) ξ | + N X i =2 σ ( e − τλ i − − e − tλ i λ i + σ " τ + N X i =2 − e − τλ i λ i . Let u i be the eigenvector associated with λ i . It follows from (5.11) that e − t B u = u , e − t B u i = e − tλ i u i → t → ∞ , i = 2 , , · · · , N, together with (5.12), which implies(5.13) ( CV ) = 1 N lim t →∞ E [ | ξ ( t + τ ) − ξ ( t ) | ] = 1 N σ " τ + N X i =2 − e − τλ i λ i . (cid:3) Concluding remarks
To model the (synchronized) beating of cardiac muscle cells, we proposed and investigated the stochasticphase equations with the irreversibility of beating (reflective boundary), induced beating and refractory.We also develop some new analysis of the conventional Kuramoto model. The application of our modelsto reproducing the bio-experimental results had been carried out in [14]. This paper mainly focuses onthe theoretical analysis, where intend to reveal the relationship between the parameters of the model andthe statistic properties of the (synchronized) beating intervals.One interesting discovery of the single-isolated cell’s model is that the distribution of beating intervalhas the coefficient variance with an upper bound p / ≈ . N -cells network.We mention some possible modifications and extensions for the proposed models, for example, thephase-dependent noise strength σ ( φ ) with σ (0) = 0, the non-interaction with other cells during refractory(i.e., A i,j = 0 for 0 ≤ φ i ≤ B i ), the irreversibility for both φ = 0 and φ = φ , and so on. Moreover, forlarge-size network, to model the propagation of the potential action (beating) of heart tissue, one canintroduce a tiny time-delay η ( η ≪
1) of the induced beating, that is if cell i beat spontaneously at time t and the neighboring cells are our of refractory, then the neighboring cells are induced to beat at time t + η . Acknowledgments
The authors would like to thank Kenji Yasuda for valuable comments. A part of this work is supportedby Core Research for Evolutional Science and Technology (CREST) of the Japan Science and TechnologyAgency (JST), Japan, and by Platform for Dynamic Approaches to Living System from the Ministry ofEducation, Culture, Sports, Science and Technology, Japan.
STOCHASTIC PHASE MODEL FOR CARDIAC MUSCLE CELLS 27
References [1] S. Abramovich-Sivan and S. Akselrod. A pacemaker cell pair model based on the phase response curve.
Biol. Cybern. ,79:77–86, 1998.[2] A. N. Burkitt. A review of the integrate-and-fire neuron model: I. homogeneous synaptic input.
Biol. Cybern. , 95:1–12,2006.[3] E. C¸ inlar.
Introduction to Stochastic Processes . Dover Publications, Inc., 2013.[4] Yu-Chuan Chang and Jonq Juang. Stable synchrony in globally coupled integrate-and-fire oscillators.
SIAM J. Appl.Dyn. Syst. , 7:1445–1476, 2008.[5] C. E. Dangerfield, D. Kay, and K. Burrage. Stochastic models and simulation of ion channel dynamics.
Procedia ComputerScience , 1(1587–1596), 2012.[6] R. L. DeHaan and R. Hirakow. Numerical simulations of angiogenesis in the cornea.
Exp. Cell Res. , 70:214–220, 1972.[7] L. C. Evans.
An Introduction to Stochastic Differential Equations . American Mathematical Society, 2013.[8] R. FithHugh. Impulses and physiological states in theoretical models of nerve membrane.
Biophysical J. , 1:445–466,1961.[9] K. Goshima and Y. Tonomura. Synchronized beating of embryonic mouse myocardial cells mediated by cells in monolayerculture.
Exp. Cell Res. , 56:387?–392, 1969.[10] M. R. Guevara and T. J. Lewis. A minimal single-channel model for the regularity of beating in the sinoatrial node.
Chaos , 5:174–183, 1995.[11] I. Harary and B. Farley. In vitro studies on single beating rat heart cells. ii. intercellular communication.
Exp. Cell Res. ,29:466–474, 1963.[12] J. M. Harrison.
Brownian Motion and Stochastic Flow Systems . John Wiley & Sons, 1985.[13] A. Hatano, J. Okada, T. Washio, T. Hisada, and S. Sugiura. A three-dimensional simulation model of cardiomyocyteintegrating excitation-contraction coupling and metabolism.
Biophys. J. , 101:2601–2610, 2011.[14] T. Hayashi, T. Tokihiro, H. Kurihara, and K. Yasuda. Community effect of cardiomyocytes in beating rhythms isdetermined by stable cells.
Scientific Reports , 7((1)):15450, 2017.[15] A. L. Hodgkin and A.F. Huxley. A quantitative description of membrane current and its application to conduction andexcitation in nerve.
J. Physiol. , 117:500–544, 1952.[16] H. J. Jongsma, L. Tsjernina, and J. deBruijne. The establishment of regular beating in populations of pacemaker heartcells. a study with tissue-cultured rat heart cells.
J. Mol. Cell Cardiol. , 15:123–133, 1983.[17] T. Kaneko, K. Kojima, and K. Yasuda. Dependence of the community effect of cultured cardiomyocytes on the cellnetwork pattern.
Biochem. Biophys. Res. Commun. , 356:494–498, 2007.[18] J. Keener and J. Sneyd.
Mathematical Physiology . Springer-Verlag, New York, 1998.[19] J. P. Keener, F. C. Hoppensteadt, and J. Rinzel. Integrate-and-fire models of nerve membrane response to oscillatoryinput.
SIAM J. Appl. Math. , 41:503–517, 1981.[20] K. Kojima, T. Kaneko, and K. Yasuda. Role of the community effect of cardiomyocyte in the in the entrainment andreestablishment of stable beating rhythms.
Biochem. Biophys. Res. Commun. , 351:209–215, 2006.[21] H. Kori, Y. Kawamura, and N. Masuda. Structure of cell networks critically determines oscillation regularity.
J. Theor.Biol. , 297:61–72, 2012.[22] Y. Kuramoto.
Chemical Oscillations, Waves, and Turbulence . Springer-Verlag, New York, 1984.[23] P. L. Lions and A. S. Sznitman. Stochastic differential equations with reflecting boundary conditions.
Comm. Pure.Appl. Math. , 37:511–537, 1984.[24] N. H. Lovell, S. L. Cloherty, B. G. Celler, and S. Dokos. A gradient model of cardiac pacemaker myocytes.
Prog. Biophys.Mol. Biol. , 85:301–323, 2004.[25] H. P. Mckean.
Stochastic Integrals . Academic Press, 1969.[26] R. M. H. Merks and P. Koolwijk. Synchronization of electrically induced calcium firings in self-assembled cardiac cells.
Biophys. Chem. , 116:33–39, 2005.[27] D. C. Michaels, E. P. Matyas, and J. Jalife. Dynamic interactions and mutual synchronization of sinoatrial node pace-maker cells a mathematical model.
Circ. Res. , 58:706–720, 1986.[28] Renato E. Mirollo and Steven H. Strogatz. Synchronization of pulse-coupled biological oscillators.
SIAM J. Appl. Math. ,50:1645–1662, 1990.[29] C. C. Mitchell and D. G. Schaeffer. A two-current model for the dynamics of cardiac membrane.
Bull. Math. Biol. ,65:767–793, 2003.[30] J. D. Murray.
Mathematical Biology . Springer-Verlag, Verlin Heiderberg, 2002. 3rd ed.[31] J. Nagumo, S. Arimoto, and A. Yoshizawa. An active pulse transmission line simulating nerve axon.
Proc. IRE , 50:2061–2072, 1962.[32] C. S. Peskin.
Mathematical Aspects of Heart Physiology . Courant Institute of Mathematical Sciences, New York Univer-sity, New York, 1975.[33] V. S. Petrov, G. V. Osipov, and J. A. K. Suykens. Influence of passive elements on the dynamics of oscillatory ensemblesof cardiac cells.
Phys. Rev. E , 79:046219:13pp, 2009. , TATSUYA HAYASHI AND TETSUJI TOKIHIRO [34] L. Sacerdote and M. T. Giraudo. Stochastic Integrable and Fire Models: A Review on Mathematical Methods and TheirApplications . Springer-Verlag, Berlin Heidelberg, 2013. in “Stochastic Biomathematical Models with Applications toNeuronal Modeling”, pp.99–148.[35] A. V. Skorokhod. Stochastic equations for diffusion processes in a bounded region.
Theory Probab. Appl. , 6:264–274,1961.[36] L. Slominski. Some remarks on approximation of solutions of SDE’s with reflecting boundary conditions.
Math. Comp.Simulat. , 38:109–117, 1995.[37] V. Torre. A theory of synchronization of heart pace-maker cell.
J. Theor. Biol. , 61:55–71, 1976.[38] A. T. Winfree.
The Geometry of Biological Time . Springer-Verlag, New York, 2001.[39] Y. Yamauchi, A. Harada, and K. Kawahara. Changes in the fluctuation of interbeat intervals in spontaneously beatingcultured cardiac myocytes: experimental and modeling studies.
Biol. Cybern. , 65:147–154, 2002. Institute of Fundamental and Frontier Sciences, University of Electronic Science and Technology ofChina, Graduate School of Information Science and Technology, Hokkaido University, The GraduateSchool of Mathematical Sciences, The University of Tokyo
E-mail address : wind [email protected], [email protected],3