Community effect of cardiomyocytes in beating rhythms is determined by stable cells

The community effect of cardiomyocytes was investigated in silico by the change in number and features of cells, as well as configurations of networks. The theoretical model was based on experimental data and accurately reproduced recently published experimental results regarding coupled cultured cardiomyocytes. We showed that the synchronised beating of two coupled cells was tuned not to the cell with a faster beating rate, but to the cell with a more stable rhythm. In a network of cardiomyocytes, a cell with low fluctuation, but not a hight frequency, became a pacemaker and stabilised the beating rhythm. Fluctuation in beating rapidly decreased with an increase in the number of cells (N), almost irrespective of the configuration of the network, and a cell comes to have natural and stable beating rhythms, even for N of approximately 10. The universality of this community effect lies in the fluctuation-dissipation theorem in statistical mechanics.

Synchronisation of biological cycles is indispensable to life 1,2 . The heartbeat is the representative phenomenon of synchronisation in physiology in which spontaneous pulsations of cardiomyocytes are tuned to a certain beating rate. Extensive research has been devoted to understanding the mechanism of regularity in beating of cardiac cells experimentally and theoretically [3][4][5][6][7][8][9] . Contraction of a cardiomyocyte is caused by complex electrophysiological processes. Detailed analyses of these processes requires elaborated mathematical models composed of a large number of equations 10,11 . However, to understand synchronisation, a small number of simultaneous ordinary equations of membrane currents and action potentials, such as the Hodgkin-Huxley equation or its reduced form, the FitzHugh-Nagumo equation and the Van der Pol equation, are sufficient for capturing the main phenomenon of cellular dynamics (see, for example 12,13 ). Most mathematical models for interacting cardiomyocytes are based on these equations [14][15][16][17][18] . A network of cardiomyocytes is regarded as a system of interacting, self-sustained (nonlinear) oscillators. To explain the details of synchronisation in such a system of oscillators, phase equations have been extensively and successfully used 19,20 . A variety of studies using phase equations in a network of cardiac cells has been reported, such as synchronisation of cardiac pacemaker cells to external periodic stimuli, phase resetting properties of cardiac cells 21 , and oscillation regularity depending on the cell networks 22 . Analysis methods of cellular neural networks have been extensively studied and some important and interesting results have been obtained [23][24][25] . Hamada et al. 26 investigated the spontaneous order in synchronisation of beating. They modelled a cell by the WJG model. This model is expressed as fairly complex, simultaneous, differential equations in terms of the membrane potential, ion currents and ion concentrations, and its extension with stochastic processes. They showed that the time interval necessary for synchronisation is strongly dependent on the strength of cell-to-cell conductance, and it is shortened by stochastic fluctuations. Fluctuation of interbeat intervals decreases as the size of the cell cluster increases, which is consistent with the experimental data.
Cells acquire a function of a tissue by forming a group. An example of this situation is that the sinoatrial node, which generates the stable beating rhythm of the normal heart, is composed of cardiomyocytes beating autonomously. This means that the function of cardiomyocytes changes from a cellular level to an organizational level by the effect of the assembly. Investigation into the process of acquiring function helps us to understand tissue models. Recently, an on-chip single-cell-based culture system was developed. Small artificial networks of cardiomyocytes can be constructed and their spontaneous beating rhythms can be measured in terms of the effects of number of cells, configurations, and types of cells 27,28 .
Although isolated cardiomyocytes are relatively heterogeneous and their beating rhythms are inconsistent, even a pair of cardiomyocytes tend to synchronise when connected with each other. Features of an individual beating cell are currently measurable and a various configuration of a cellular network, which greatly affects cell-to-cell interactions, is artificially constructible. Therefore, examining how heterogeneity of cells and cell-to-cell interactions affect synchronisation in a small cluster of cardiomyocytes would be of considerable importance. There are two important observable quantities in a cardiomyocyte; one is its cell cycle (beating rate) and the other is its refractory period. In particular, a cardiomyocyte has a relatively long refractory period compared with that of a neuron. The difference in refractory periods among cardiomyocytes is expected to affect the behavior of synchronisation.
When cardiomyocytes are isolated, they only beat independently. However, if cardiomyocytes come into contact and interact with each other, their beating rhythms become synchronised. Researchers originally hypothesised that, in a network of cardiomyocytes, firing of one cardiomyocyte triggers induced firing of adjacent cardiomyocytes, and all of the cardiomyocytes start beating synchronously, and that the beating rate is tuned to the fastest cardiomyocyte 5 . However, recent experiments have shown that other cells are synchronised not to the fastest cell, but to the cell with the least fluctuation in beating rhythm 29 (Supplementary Fig. S1).
The present study aimed is to investigate the community effect of cardiomyocytes in different configurations of networks constituted by cells with specified characteristics of beating rhythms. We also aimed to clarify how an assembly of cells acquires stability, one of the most important universal features in biological systems. Because preparing a cardiomyocyte with given properties in an in vitro experiment is relatively difficult, we developed a mathematical model, which explained this behavior of cardiomyocytes with high reliability. The model that we adopted in this study is a modification of the well-known integrate-and-fire model, which has been widely used as a spiking neuron model 30,31 . Using this model, we investigated the dependence of fluctuation in beating rhythm on the number of cells and that on configuration of the networks.

Methods
We constructed a theoretical modelling for synchronisation of cardiomyocytes by modifying the integrate-and-fire model. The integrate-and-fire model disregards biological details and only focusses on causal relationships in the phenomenon. This model consists of biological oscillators or, equivalently, phase variables just as those in phase equations, but the interaction between them is instantaneous and spiky. Based on the simple Peskin's model 32 , we included refractory periods, stochastic process, and weak cell-to-cell interactions, which modulate phase variables 19,20 . An important point is that the stochastic process and the cell-to-cell interaction are correlated through the fluctuation-dissipation theorem that gives the relation between fluctuations and linear response to external force 33,34 . Despite the simplicity of our model, it accurately reproduces most experimental results for interaction between two cardiomyocytes, though it has only one free parameter and the other parameters are determined by experiments. In the current study, we used data for 14 pairs of cardiomyocytes reported from a previous publication 29 . In the experiment 29 , cardiomyocytes that were dissociated from 13 to 14-day-old mouse embryos were used. We used Fortran and Mathematica to perform numerical simulations and to analyze the obtained data.

Mathematical modelling.
When considering a network of N cardiomyocytes, the model was described by , which denote the state of ith cardiomyocytes. We assumed that the ith cardiomyocyte fires (beats) when φ π = ≡ t ( ) 0 ( 2 ) i . This firing occurs either at φ t ( ) i , reaches 2π, or the following conditions are satisfied: Additionally, one of the cardiomyocytes connected to ith cardiomyocyte (e.g., jth cardiomyocyte) fired at a delayed time τ ago (i.e., ). Otherwise, we assumed that t ( ) i φ is governed by the following interacting stochastic differential equation. Our mathematical modelling for cell-i is as follows: where ω i is the average phase velocity of cell-i, θ i is a phase corresponding to the refractory period of cell- is a stochastic process with standard deviation σ, and Σ j denotes the summation over all the cardiomyocytes connected with the ith cardiomyocyte. φ φ V( , ) i j denotes the weak interaction through the membrane potential, which we assumed as the following form: Here μ is a positive constant. Note that ω i , θ i , σ i can be determined by single-cell experiments for each cardiomyocyte. We used an extended random walk as the stochastic process W( ) σ . The positive constant μ is the only parameter in our model that cannot be directly measured by experiments. The dynamics of the state variable t ( ) i φ is schematically shown in Fig. 1.
Numerical simulation method. In our simulation, the stochastic process is expressed by an extended random walk and the standard deviation is defined by σ = Δ Δ x t : / 2 , where Δt and Δx are the time interval and the one-step distance in the random walk, respectively. We used the following difference equations as a numerical approximation of equation (1). For almost all cardiomyocytes, we considered an ordinary random walk as follows: where Δt is the time difference interval for numerical simulations, the retardation time τ is set as Δt × m nonnegative integer, and x t 2 is the spatial difference determined from σ i . However, for cardiomyocytes with a large fluctuation, we could not reproduce the same fluctuation in beating rhythm by an ordinary random walk. This is because the coefficient variation (CV%), which was defined by 100× standard deviation/mean beating rate, could be less than  . 100 2/3 81 65 35 . However, some cardiomyocytes have the CV% which exceed this value. Therefore, we adopted the following extended random walk, which is a history-dependent stochastic process, when beating fluctuation was larger than 81.65: However, By choosing appropriate values of q and r, we could reproduce the large fluctuation observed in the experiments.

Results
Comparison of the model with experimental results of two cardiomyocytes. In the previous experiments using cultured cardiomyocytes 29 , the mean beating rate and its fluctuation before and after synchronisation were observed for 14 pairs of cardiomyocytes. We applied our mathematical model of the present study to determine whether it could reproduce the results of these pairs of cardiomyocytes. We numbered these 14 pairs from Nos1 to 14 and distinguished the two cardiomyocytes in a pair by denoting "cell-1" and "cell-2". For each  (1) for cell-i (i = 1, 2), so that the model reproduced the same mean beating rate and fluctuation in beating rhythm. Refractory periods of cardiomyocytes are almost the same as those for normal cells. Therefore, we assumed that each cell had the common refractory period The exact values of parameters , , are shown in data tables in the Supplementary Information (Supplementary Table S1). Figure 2 shows the mean beating rates and fluctuation in the beating rhythm after synchronisation for the 14 pairs obtained by the experiments and numerical simulation by our model. We used the retardation time τ = 0 because it was estimated as 10 −3~1 0 −4 of the mean beating rate. We used μ = 6.5 in numerical simulations. By defining an evaluation function, we were able to determine the free parameter μ to minimize the function. The results of the numerical simulations were almost constants for a relatively wide range of μ. This finding indicated that our model was robust against the free parameter μ. The dependence of theoretical calculation on μ is shown in the Supplementary Information (Supplementary Note S1 and Supplementary Fig. S2). Except for pair No. 14, the simulated values accurately agree with the experimental values. Fluctuation in beating of a pair of synchronised cardiomyocytes almost coincided with that of less fluctuating cardiomyocytes, while the mean beating rate after synchronisation was widely distributed. Some synchronised cardiomyocytes coincided with faster rates, some with slower rates, and others with intermediate rates. The experimental result of pair No.14 is exceptional because it is the only pair in which fluctuation increased after synchronisation.

Comparison with the Kuramoto model. For two coupled oscillators
2 . For simplicity, we assumed that f(x) is 1-periodic function with f(−x) = −f(x) (e.x. f(x) = sin(x)) and A i,j ≥ 0. The two-oscillators phase model (Kuramoto model) with noise is as follows: for i, j = 1, 2, i ≠ j, where ω i and σ i denote the drift and noise strength constants, respectively, and = W { } i i 1,2 is independent standard Brownian motion. For two cases (Case (i) and Case (ii)), we applied the Kuramoto model (8) and our model (1) to synchronisation of two coupled cardiomyocytes. The numerical simulation results were compared with biological experiment data 29 .
Case (i) A case of synchronisation to a cardiomyocyte with a fast and stable beating rhythm. When cell-1 which had a mean beating rhythm of 0.64 s and fluctuation of 12.3 [CV%] and cell-2 with 1.23 s and 25.1 [CV%] were coupled, we found that the fast and stable cardiomyocyte (cell-1) acted as a pacemaker and the beating rhythm after synchronisation was tuned to cell-1 (Fig. 3a). We investigated whether our model and the Kuramoto model could reproduce the experimental results. Figure 3b and c shows the theoretical predictions from our model and the Kuramoto model, respectively, and both models reproduced the observations. The mean beating rate and beating fluctuation for the experimental result, our model, and the Kuramoto model are shown in the Supplementary Information (Supplementary Table S3).

Case (ii) A case of synchronisation to a cardiomyocyte with a slow and stable beating rhythm.
When cell-1 which had a mean beating rhythm of 1.10 s and fluctuation of 149 [CV%] and cell-2 with 1.40 s and 41.2 [CV%] were coupled, we found that the slow and stable cardiomyocyte (cell-2) acted as a pacemaker and the beating rhythm after synchronisation was tuned to cell-2 (Fig. 3d). When we compared the numerical result of our model with  Table S2).
that of the Kuramoto model, our model was closer to the experimental data than the Kuramoto model. Our model showed that the beating rhythm after synchronisation was tuned to the slow and stable cardiomyocyte (Fig. 3e). However, the Kuramoto model showed that beating fluctuation of the slow and stable cardiomyocyte was increased after synchronisation, which differed from the observations (Fig. 3f). The mean beating rate and beating fluctuation of the experimental result, those of our model and those of the Kuramoto model are shown in the Supplementary Information (Supplementary Table S4). Therefore, our model showed that even though the mean beating rate of a cardiomyocyte was slow, a cardiomyocyte with more stable beating fluctuation dominated the beating rhythm after synchronisation. If CV values were set to a similar values, we can also obtain the same results (Supplementary Note S2, Supplementary Fig. S3 and Supplementary Table S9). In previous numerical simulations, we did not consider the effect of retardation  time (τ = 0). When we incorporated this effect, the behavior of our model barely changed because of the existence of a refractory period much longer than τ. However, if the refractory period is not taken into account, then a couple of cardiomyocytes continuously fire with the period of the retardation time, which is biologically unacceptable. In a system of two cardiomyocytes, we can use a retardation time τ = 0, but we should consider the effect of retardation time as the system size increases. In this case, the existence of the refractory period will have significant effects on the system.
As an application of our mathematical modelling, we then performed two numerical experiments on networks of cardiomyocytes and investigated the community effect of cardiomyocytes. = .
= . = .  Size and configuration dependence on fluctuation of the system. First, we investigated the dependence of fluctuation in beating rhythm of cardiomyocytes on the size and configuration of the system. The configurations that we considered were star, 2D lattice and 1D lattice networks (Fig. 4). . In all configurations, fluctuation rapidly decreased with an increase in the size of the system. Among the three configurations, a reduction in fluctuation was most rapid in the 2D lattice network, and fluctuation in the 1D lattice network was always larger than that in the other two configurations. Furthermore, we considered the larger size (about 1000 cells) of the network in the 2D lattice network. Figure 5d-f shows the size dependence of fluctuation of the 2D lattice network. Similar to Fig. 5a-c, all the elements had the same beating properties. In all cases, the beating fluctuation decreased as the community size increased. In Fig. 5(a-f), beating fluctuation of model cardiomyocytes are much different, but if CV values were set to a similar values, the same results was obtained ( Supplementary  Figs S4 and S5). For an ordinary stochastic ensemble, such as an independently identical distributed ensemble, the dependence of standard deviation of fluctuation on system size N was proportional to N −1/2 . However, the data of fluctuation plotted on the graph (Fig. 5) considerably deviated from the line of N −1/2 and the feature of beating fluctuation was relatively different from that of ordinary stochastic ensembles.
Dependence of cellular properties and numbers on fluctuation of the system. We then investigated the change in beating rhythms after connecting two subsystems of cardiomyocytes. First we composed  . Schematic explanation of why the beating rhythm tend to be synchronised to that of more stable cardiomyocytes after connection of two cardiomyocytes. A stable cardiomyocyte can be compared with a heavy pendulum and an unstable cardiomyocyte with a light pendulum. (a) External fluctuation has little effect on a pendulum's period of swing if its weight is heavy, but has strong effects if its weight is light. (b) When two pendulums are coupled and synchronised, their period of swing is close to that of the heavier pendulum, and fluctuation will be reduced because the total mass of weight increases. referential subsystems of four model cells and nine model cells so that these subsystems had the property of a standard beating rhythm (mean beating rate 1.20~1.30 s; fluctuation 15.0~20.0 [CV%]) As for the subsystems which are connected to referential subsystems, we considered subsystems consisting of four types of cardiomyocytes: (i) first and stable cell, (ii) first and unstable cell, (iii) slow and stable cell, and (iv) slow and unstable cell. For networks, we considered the three types of configurations shown in Fig. 4. A single cell was connected to a centre cell of the referential star network, to a cell on a link of the referential 2D lattice network, and to a cell on an edge of the 1D lattice network (Fig. 6a,b). In case of subsystems of the same size, we connected them by a single link between two cells at the same positions in the configuration. We used the two centre cells for the star networks, cells on the links for the 2D lattice network, and the two cells at the edges for the 1D lattice network (Fig. 6c,d).
The results are shown in the Supplementary Information (Supplementary Tables S5-S8). The mean beating rate and fluctuation (CV%) of each cell were estimated by the data of 3,000 firings. We showed three typical results of the numerical simulation. First, we considered the referential 2D lattice network with nine cells and a single cell that had a fast and unstable beating rhythm (Fig. 7a). In this case, the beating rhythm of the system after connection was tuned to the rhythm of nine cells with more stable beating (Fig. 7b). Second, the rhythm of a single cell was fast and stable (Fig. 7c). We then found that even a single cell could lower fluctuation of the referential network and the synchronised beating rhythm was tuned to a stable single cell (Fig. 7d). Finally, we considered a referential 2D lattice network with four cells and four cells grouped with a fast and stable beating rhythm (Fig. 7e). The beating rhythm after connection was also tuned to the rhythm of the more stable group (Fig. 7f). In the above three cases (Fig. 7a-f), every cell started synchronising after connection and fluctuation of the cells became equal in the combined system. However, in the case of 1D lattice network in Fig. 6d, the cells near the edge of 1D lattice network with nine cells showed an exceptionally large fluctuation compared with the other cells. Synchronisation did not occur in the referential 1D lattice network with nine cells + nine cells with a fast and stable beating rhythm. Furthermore, fluctuations of combined systems reduced their intensity, except for when there was a single cell or cell group with a slow and stable property. However, the increment in fluctuation was small, even in these cases (e.g., Fig. 7g,h). When a referential subsystem was connected to a counterpart consisting of one of the other three types of cells, the constituent cells acquired a common intensity of fluctuation. The intensity was intermediate between that of the prior two subsystems, but was similar to that of the less fluctuating subsystem.

Discussion
We investigated the community effect of networks of cardiomyocytes by using an interacting integrate-and-fire model with a refractory period. The reliability of the present model was verified by accurately reproducing recent experiments on pairs of cultured cardiomyocytes by Kojima et al. 29 , despite the fact that the model has only one tuning parameter. Interestingly, the fluctuation observed in their experiments cannot be explained by simple Brownian motion or equivalently random walks. This is because some of the fluctuations in beating rhythm (CV%) exceeded the theoretical limit evaluated for Brownian motion. An important observation in their experiments is the finding that a pair of cardiomyocytes, when connected, tended to beat synchronously at a rate of the cell with a stable beating rhythm, but not the cell with a faster beating rhythm. This community effect of cardiomyocytes towards stability was confirmed with the present model by examining networks of various configurations and constituent cells with various beating rhythms. Even a single stable cardiomyocyte could lower fluctuation of a network consisting of some cardiomyocytes. The reason why an unstable cardiomyocyte tends to follow a stable cardiomyocyte may be explained as follows. A stable cardiomyocyte has the property where its dynamics are only slightly affected by external or internal disturbance. Therefore, there is little effect of interactions from neighboring cardiomyocytes, while an unstable cardiomyocyte has the opposite property, and is strongly affected by its neighbors. A stable cardiomyocyte corresponds to a pendulum with a heavy mass in contrast to an unstable cardiomyocyte that corresponds to that with a light mass (Fig. 8). When we connect these pendulums, clearly the pendulum with a light mass tends to follow that with a heavy mass. This feature is a consequence of the fluctuation-dissipation theorem, which provides a universal relation between fluctuation and a linear response 33,36 . In our model, the factor i 2 σ of the interaction term φ φ ∑ V( , ) j j i in equation (1) was due to this theorem, and it plays an essential role in stabilising the beating rhythm after synchronisation. Stability is one of the most important and universal features of biological systems. Interestingly, one of the origins of biological stability is a universal principle in statistical physics, that is, fluctuation-dissipation theorem.