Delayed global feedback in the genesis and stability of spatiotemporal patterns in paced biological excitable media

A multi-scale approach was used to investigate the roles of delayed global feedback (DGF) in the genesis and stability of spatiotemporal patterns in periodically-paced excitable media. Patterns that are temporal period-2 (P2) and spatially concordant (in-phase) or discordant (out-of-phase) were investigated. First, simulations were carried out using a generic spatiotemporal model composed of coupled FitzHugh-Nagumo units with DGF. When DGF is absent, concordant and discordant P2 patterns occur depending on initial conditions. The discordant P2 patterns are spatially random. When the DGF is negative, only concordant P2 patterns exist. When the DGF is positive, both concordant and discordant P2 patterns can occur. The discordant P2 patterns are still spatially random, but they satisfy that the global signal exhibits a temporal period-1 behavior. Second, to validate the spatiotemporal dynamics in a biological system, simulations were carried out using a 3-dimensional physiologically detailed ventricular myocyte model. This model can well capture the intracellular calcium release patterns widely observed in experiments. The properties of DGF were altered by changing ionic currents or clamping voltage. The spatiotemporal pattern dynamics of calcium release in this model match precisely with those of the generic model. Finally, theoretical analyses were carried out using a coupled map lattice model with DGF, which reveals the instabilities and bifurcations leading to the spatiotemporal dynamics and provides a general mechanistic understanding of the role of DGF in the genesis, selection, and stability of spatiotemporal patterns in paced excitable media. Author Summary Understanding the mechanisms of pattern formation in biological systems is of great importance. Here we investigate the dynamical mechanisms by which delayed global feedback affects pattern formation and stability in periodically-paced biological excitable media, such as cardiac or neural cells and tissue. We focus on the formation and stability of the temporal period-2 and spatially in-phase and out-of-phase patterns. Using a multi-scale modeling approach, we show that when the delayed global feedback is negative, only the spatially in-phase patterns are stable; when the feedback is positive, both spatially in-phase and out-of-phase patterns are stable. Also, under the positive feedback, the out-of-phase patterns are spatially random but satisfy that the global signals are temporal period-1 solutions.


Introduction
DGF in the genesis of spatiotemporal dynamics in an array of coupled FHN units 136 We used a generic model consisting of a one-dimensional (1D) array of coupled FHN units 137 to investigate the spatiotemporal excitation patterns. The governing differential equations are: In Eq.3, n is the index of the beat number, T is the pacing period, T is the pulse duration, and  147 is the feedback strength. is the peak value of the spatial average of (denoted as ) at the

Excitation patterns without DGF
In the absence of DGF, i.e., in Eq.3, a bifurcation from temporal P1 to P2 occurs as = 0 154 the pacing period T decreases (Fig.S2). We found that when T first passes through the bifurcation 155 point, the system can only exhibit a spatially concordant P2 (Con-P2) pattern; as T decreases 156 further, the system can exhibit a Con-P2 (Fig.2a) or a spatially discordant P2 (Dis-P2) pattern 157 (Fig.2b), depending on initial conditions. It appears that the probability of forming a Dis-P2 pattern 158 increases as the spatial heterogeneity of the initial condition increases (Fig.2c). Moreover, the Dis-159 P2 patterns are spatially random and selected by initial conditions. To quantify this property, we 160 measured the spatial domain sizes (see Fig.2b for definition) from 2000 random trials for a given 161 standard deviation of the spatial heterogeneity of the initial condition, and plotted the 162 corresponding histogram (Fig.2d). It shows that the domain size can be any value as long as it is 163 greater than a minimum domain size l min , i.e., the domain sizes distribute between l min and L-l min .

164
Because of this randomness in pattern selection, the corresponding histogram of the global P2 165 amplitude ( as defined in Fig.2a) also exhibits a continuous distribution (Fig.2e). Δ 168 To investigate the effects of DGF on the spatiotemporal pattern dynamics, we carried out 169 simulations by scanning the pacing period T and DGF strength  (Fig.3a). There are four distinct 170 regions: uniform P1 pattern (yellow), Con-P2 pattern only (cyan), Dis-P2 pattern only (black), and 171 both concordant and discordant P2 (Con/Dis-P2) patterns (red). The blue curve is the stability 172 boundary between P1 and P2 for a single uncoupled FHN unit. For , only uniform P1 and < 0 173 Con-P2 patterns were observed, independent of initial conditions. The uniform P1 and Con-P2 indicating that the dynamics in the 1D array is the same as in the single FHN unit. For , a > 0 transition from uniform P1 to Dis-P2 occurs as T decreases (from yellow to black), which is caused 177 by a spatial-mode instability of the uniform P1 state. As T decreases further (red region), both 178 Con-P2 and Dis-P2 patterns can occur depending on the initial conditions (Fig.3b). 179 Furthermore, we performed the same statistical analysis as in the case of no DGF (Fig.2 d   180 and e) for different regions. In the Dis-P2 only region (Fig.3c), the domain sizes distribute between 181 0 to L/2 (more accurately, the domain size can be L/2 and any value between l min and L/2-l min ), but 182 remains zero for all patterns. In the Con/Dis-P2 region (Fig.3d), the distributions are similar model undergoes a bifurcation from P1 to P2 (alternans) as the pacing period T decreases (Fig.S2).

202
We investigated the subcellular Ca 2+ release patterns under both AP clamp and free-running  Ca 2+ -to-APD coupling properties on Ca 2+ release patterns and then linked them to DGF. 227 We systematically explored the spatiotemporal dynamics by altering the pacing period T 228 and the maximum conductance of the two currents, as summarized in Fig.5a. When the Ca 2+ -to-229 APD coupling is negative (large I SK ), a transition from uniform P1 to Con-P2 patterns occurs as T 230 decreases, and this transition occurs at a larger T value as the maximum I SK conductance increases.

231
When the coupling is positive (large I nsCa ), a transition from uniform P1 to Dis-P2 patterns (yellow 232 to black) occurs as T decreases. Under both coupling conditions, as T decreases further, the system 233 enters the Con/Dis-P2 regime (red), in which both Con-P2 and Dis-P2 patterns can occur 234 depending on initial conditions (Fig.5b). However, as T decreases even further, the Con/Dis-P2 235 regime switches into a Dis-P2 only regime when the Ca 2+ -to-APD coupling is negative (large I SK ) 236 and into a Con-P2 only region when the Ca 2+ -to-APD coupling is positive (large I nsCa ). Therefore, 237 for the same Ca 2+ -to-APD coupling, as T decreases, the spatiotemporal patterns change from Con-238 P2 only to Dis-P2 only through a Con/Dis-P2 region or in reverse order depending on the coupling 239 properties.

240
To reveal the statistical properties of the Dis-P2 patterns, we show the histograms of 241 domain sizes and the whole-cell alternans amplitude for a parameter point in the Dis-P2 region 242 (Fig.5c) and a point in the Con/Dis-P2 region (Fig.5d). The domain size distributions for Dis-P2 243 patterns are continuous and the whole-cell alternans amplitudes always remain zero, indicating that the patterns are spatially random but always satisfying that the global signals are P1 solutions. 245 These behaviors are the same as in the model of a coupled array of FHN units (Fig.3). 246 To link the spatiotemporal Ca 2+ dynamics to the DGF properties, we performed an analysis 247 to reveal the DGF properties and their relationship with the Ca 2+ -to-APD coupling properties. The  mode is the same as that of a single uncoupled unit. Since for in Eq.7 does not depend on > 0 292 , then the feedback has no effects on the stability of the uniform P1 state for non-zero mode.

293
Because of this, the stability boundary separating uniform P1 from Dis-P2 appears to be a 294 horizontal line independent of  (Fig.6a, solid).   Fig.6a), indicating that stable Dis-P2 patterns can only exist when . Histograms of domain > 0 315 size and example spatial patterns from three locations marked in Fig.6a are plotted (Fig.6 b and c).      Histogram of global P2 amplitude from the same simulations. was measured using ∆ Δ 653 the last two beats. d. Same as panel c but T=350 ms and . = 3.5