## Abstract

The muscular layer of the uterus (myometrium) undergoes profound changes in global excitability prior to parturition. Here, a mathematical model of the myocyte network is developed to investigate the hypothesis that spatial heterogeneity is essential to the transition from local to global excitation which the myometrium undergoes just prior to birth. Each myometrial smooth muscle cell is represented by an element with FitzHugh–Nagumo dynamics. The cells are coupled through resistors that represent gap junctions. Spatial heterogeneity is introduced by means of stochastic variation in coupling strengths, with parameters derived from physiological data. Numerical simulations indicate that even modest increases in the heterogeneity of the system can amplify the ability of locally applied stimuli to elicit global excitation. Moreover, in networks driven by a pacemaker cell, global oscillations of excitation are impeded in fully connected and strongly coupled networks. The ability of a locally stimulated cell or pacemaker cell to excite the network is shown to be strongly dependent on the local spatial correlation structure of the couplings. In summary, spatial heterogeneity is a key factor in enhancing and modulating global excitability.

## 1. Introduction

The myometrium is the muscular layer that constitutes the bulk of the uterine wall. It is a syncytium of interconnected smooth muscle cells, forming an excitable medium [1], i.e. a nonlinear dynamical system that can propagate signals over long distances without damping. In the myometrium, these propagating signals trigger phasic contractions [2]. The behaviour of the myometrium as an excitable medium is thought to be influenced by the spatial variations in the excitability of individual muscle cells as well as the strength of their interconnections [3]. The modulation of global excitation by network heterogeneity may play an important role in the myometrium during pregnancy. Sufficient coupling is needed for excitation to spread, but this coupling need not be uniform, or even exist between all cells. The aim of this paper is to examine how spatial heterogeneity, in particular local variations in cell connectivity, affects global excitability. Furthermore, we study the ability of pacemaker cells to drive the network as a function of its spatial heterogeneity.

In the rodent (as in most mammals), the myometrium consists of an inner circular layer and an outer longitudinal layer of smooth muscle cells [4]. In humans, these layers are less distinct [5]. Through most of pregnancy, the myometrium remains in a predominantly quiescent state as the fetus develops [6]. However, in the days leading up to parturition, contractile activity in the myometrium undergoes major changes which prepare the uterus for labour [6]. This activation phase involves molecular changes that lead to an increase in contraction frequency, when compared with mid-gestation [7]. A key characteristic of this phase is a profound change in the connectivity between the myometrial smooth muscle cells.

If a stimulus is applied to a smooth muscle cell, its membrane potential undergoes a depolarization before eventually returning to the quiescent resting value. An excitable cell exhibits an all-or-nothing response: a cell either responds with a full excursion of the membrane potential, or barely at all. In particular, when the stimulus exceeds an excitability threshold, the response is an action potential whose magnitude is independent of the size of the (suprathreshold) stimulus. A consequence of this all-or-nothing behaviour is that a stimulus of sufficient amplitude can be reliably propagated across the network.

Pacemaker cells, by contrast, do not require an external stimulus, but exhibit periodic excitations that drive neighbouring non-pacemaker cells. The presence of specialized pacemaker cells in myometrial tissue is much disputed [8–11]. If present, it is likely that the pacemaker cells are not anatomically distinct from non-pacemaker cells and have no fixed location in the uterus [12]. Rather, they would occur dispersed within the matrix of non-pacemaker cells. Therefore, the cells can only be recognized through patterns of electrical activity: a slow depolarization of the smooth muscle membrane which results in the generation of an action potential [13]. The oscillation frequency of a pacemaker cell dictates the rate of uterine contractions.

In order to generate action potentials, a myometrial cell maintains transmembrane gradients of several ionic species by means of active transport across the cell membrane. In humans, the action potential is initiated by an inward, depolarizing current carried by calcium ions [2]. Post-excitation, the cells enter a refractory period, during which they are temporarily unable to become excited [14]. Electrical coupling between myometrial smooth muscle cells is maintained by intercellular channels through which ions and certain metabolites can pass from one cell to another. Each channel consists of two connexons, one contributed by each of the communicating cells; furthermore, each connexon is a complex of six connexin proteins [15]. A gap junction is a cluster of such channels joining two cells electrically. Myometrial gap junctions vary in number over the course of gestation; in fact, an increase in the number of gap junctions in myometrial tissue is indicative of the onset of labour [1]. Moreover, the precise nature of the conditions that stimulate the expression of gap junctions is of considerable clinical importance and an understanding of these conditions may ultimately aid early diagnosis and improve management of pre-term labour [16].

Heterogeneity has long been known to be a factor in the spatial patterns of activation in excitable systems [17]. In particular, electrical heterogeneities can play a role in cardiac arrhythmias and cause a decrease in propagation velocity through the tissue. The propagation of excitation in spatially extended systems through spiral waveforms has been well studied [18,19]. Here, we focus on a question which has received much less attention: whether spatial heterogeneity in connectivity is sufficient to modulate global excitability. We use a mathematical model to investigate how the global excitation of the myometrium is affected by the spatial heterogeneity of the system. Benson *et al.* [20] analysed a heterogeneously coupled model which is based on FitzHugh–Nagumo dynamics, but this was a continuum-model in which the coupling between cells was not represented explicitly. Spatial heterogeneity was found to assist the transition between quiescence and excitability. However, heterogeneity in coupling alone was not observed to produce synchronous activity. Following these authors, we use the FitzHugh–Nagumo model [21,22] and introduce heterogeneity by generating stochastic coupling structures based on several statistical distributions. These include empirical distributions obtained from experiments performed on myometrial cells taken from mice at the 15th and 18th day of pregnancy (i.e. towards the end of gestation). In addition, we investigate how spatial heterogeneity of the intercellular connections modulates the ability of a pacemaker cell to drive the network.

## 2. Material and methods

Both *in silico* and *in vitro* methods were employed. Myometrial tissue was simulated in a mathematical network model to determine the effect of spatial variations of coupling strengths on the spread of excitation. Since coupling depends on cell capacitance, as explained later, cell capacitance values were taken from pregnant mouse myometrium to define a statistical distribution for coupling values, which was used to construct an asymmetrically coupled model. Furthermore, data from pregnant mouse myometrial cells were used to simulate statistical variation in resting membrane potential for each cell. These physiologically realistic spatial statistics were also used to determine the efficacy of pacemaker cells.

### 2.1. Experimental methods

#### 2.1.1. Electrophysiological measurements

Animals used were C57BL/6 mice which were time-mated (within a 2 h period) to generate pregnancies with an accurate gestational age. Mice were sacrificed by carbon dioxide inhalation at gestation day 15 or 18, time points towards the end of pregnancy. Strips of myometrium from the longitudinal layer (2 × 2 × 20 mm) were dissected from freshly isolated uteri in ice-cold physiological saline. The strips were washed in Ca^{2+} and Mg^{2+}-free HBSS (Fisher Scientific) at 37°C for 10, 20 and 30 min sequentially. This was followed by a 45-min incubation in HBSS containing Liberase TM (Roche) at a final concentration of 0.13 units ml^{−1}. Digestion was terminated by several dilutions with fresh HBSS. Cells were dispersed by slow trituration through a wide-bore fire polished glass pipette in HBSS solution. Single myometrial cells were filtered through a 200 μm gauze and stored in HBSS for use within 6 h.

The cell membrane was perforated using the antibiotic amphotericin B (600 μg ml^{−1}) [23]. Cell capacitance was measured in the enzymatically isolated smooth muscle cells using the membrane test facility of the Axopatch 700B amplifier (Axon Instruments). Repetitive square voltage pulses of *Δ**V _{m}* = 10 mV were applied from the holding potential of –60 mV and the current response was measured. The integral of the current as a function of time over the decay phase equals the unloaded charge

*Q*

_{c}which is related to the step change in voltage by

*Q*

_{c}=

*C*

_{m}*Δ*

*V*, whence

_{m}*C*can be calculated. Electronic supplementary material, table S1 lists the cell capacitance data obtained in our experiments.

_{m}Transmembrane potentials were recorded with an amplifier (Axopatch 700B; Axon Instruments) and a Digidata 1440a computer interface running pCLAMP 10.2 software (Molecular Devices, Sunnyvale, CA, USA). Resting membrane potential values were taken as the mean potential (mV) for the 5-s period immediately after an action potential has occurred.

#### 2.1.2. Data analysis

For the 18-day cell capacitance data, a distribution was fitted as follows. The coefficient of variation was calculated as *ξ*_{1} = *σ*_{18}/*μ*_{18}, where *μ* and *σ* are the mean and s.d. of the experimental data, respectively. The data were normalized by the average value, giving a normalized mean of 1 and a normalized s.d. of *ξ*_{1}. A lognormal distribution can be characterized by two parameters *λ* and *τ* such that the mean is equal to exp(*λ* + *τ*^{2}/2) and the variance equals exp{2*λ* + *τ*^{2}}(exp{*τ*^{2}−1}). These parameters satisfy the following system:
2.1The non-dimensional cell capacitance for 18-day pregnant mouse data was found to be distributed as e* ^{Z}*, where

*Z*follows a normal distribution with mean –0.03551 and s.d. 0.2665.

The 15-day data were analysed in a similar fashion. The coefficient of variation was defined as *ξ*_{2} = *σ*_{15}/*μ*_{15} using the mean and s.d. of the experimental data. However, the data were normalized by the 18-day mean *μ*_{18}, with normalized mean *α* = *μ*_{15}/*μ*_{18}. Hence
2.2The non-dimensional cell capacitance for 15-day pregnant mouse data was found to be distributed as e* ^{Z}*, where

*Z*follows a normal distribution with mean –0.2114 and s.d. 0.3073.

### 2.2. Myometrial network model

The myometrial network model consists of excitable elements interconnected by resistors representing gap junctions. Heterogeneity is introduced in the form of random variations in the strength of the cell-to-cell couplings. We use various statistical models to generate the networks, including two based on cell population statistics taken from the mice data.

#### 2.2.1. Cell model

The Hodgkin–Huxley model of electrogenic cell activity, proposed in 1952, forms the basis for the study of excitable systems [14,24–27]. While the Hodgkin–Huxley model has four state variables, FitzHugh pointed out that the essential dynamical properties are captured by a two-dimensional simplified model [21], for which Nagumo proposed an electric analogue circuit [22]. The simplified model is an excitation–relaxation oscillator with a fast (‘excitation’) state variable *v* that corresponds to the cell's membrane potential and a slow (‘recovery’) variable *w* that corresponds to gating kinetics which repolarize the excited cell. The following ordinary differential equations describe this two-variable model:
2.3aand
2.3bwhere *I* is the input current and *A*, *α*, *γ*, *w*_{0}, *v*_{0} and *ε* are positive parameters. The values shown in table 1 were used in the simulations presented in §§3.1–3.7, with the exception of *ε* which was scaled according to cell capacitance in §3.6.

The equations were solved using the NDSolve function in Mathematica, which uses a Livermore Solver for Ordinary Differential Equations approach. Results were numerically stable under variation in step size.

The behaviour of an isolated cell with the input *I* = 0 is shown in a phase-space diagram (figure 1*a*). Whereas the null isocline of *w* is a straight line, the null isocline of *v* is a cubic polynomial with an unstable branch in the middle, indicated as a dashed line. In the region left of this unstable branch, the phase point tends towards the left stable branch and ultimately towards the intersection of the null isoclines, which forms a stable stationary point. To the right of the unstable branch, the phase point tends rapidly towards the right stable branch, where the slow dynamics of *w* will drive it upwards until the branch point is attained. The rapid dynamics of *v* subsequently drives the phase point back to the left stable branch.

Excitation corresponds to an excursion along the right stable branch. To reach this branch from the stationary point, a perturbation (*Δ**v*) has to be applied on *v*. If this perturbation is insufficient to move the phase point beyond the unstable branch, the phase point rapidly relaxes back to the equilibrium (figure 1*b*), whereas sufficiently large perturbations trigger a substantial response (figure 1*c*). The unstable branch of the null isocline thus represents a threshold for excitation.

#### 2.2.2. Lattice model

The individual excitable cells were coupled together through resistors into an *n* × *n* lattice, as illustrated in figure 2*a*. The lattice was modelled with closed boundary conditions. The resistors represent the gap junctions between any two adjacent cells. The equivalent resistance of the gap junction can be calculated on the basis of the properties of the individual connexon channels. Let the individual connexon have a resistance *r _{i}* (figure 2

*b*). Suppose that the gap junction between a given pair of cells consists of

*n*connexons. Because the connexons conduct in parallel, Kirchoff's law for parallel resistors determines the total resistance

*R*of the gap junction (figure 2

*c*) 2.4If the connexons have equivalent resistance,

*r*≡

_{i}*r*and

*R*=

*r*/

*n*.

#### 2.2.3. Coupling constants

A coupling constant *K* is assigned to each gap junction and defined as follows. The cell membrane is modelled as a capacitor in parallel with a resistor, as indicated in figure 2*d*. Let *Q _{i}* denote the membrane charge of cell

*i*,

*C*the cell capacitance and

_{i}*V*the membrane potential. These quantities are related by 2.5The gap junctional current between two cells

_{i}*i*and

*j*is given by (

*V*–

_{j}*V*)/

_{i}*R*. Define the coupling constant between cells

_{ij}*i*and

*j*as

*K*= (

_{ij}*C*)

_{i}R_{ij}^{−1}and consider a cell connected to four other cells in a rectangular grid, as shown in figure 2

*d*. The gap junctional current for the central cell is a sum of four gap junctional currents where (

*i*,

*j*) denotes the location of the cell on the grid. Hence the voltage dynamics for cell (

*i*,

*j*) is given by 2.6where

*I*

_{Ch}is the total current carried by the ion channels of cell (

*i*,

*j*). The coupling values between any two adjacent cells are as follows: 2.7The coupling value is rendered dimensionless in accordance with the scaled equations for dynamics of the cell network (equations (2.3

*a*,

*b*). As defined here,

*K*is a rate constant. Dimensionless time

*t*in the FitzHugh–Nagumo model is defined by

*t*=

*R*

_{1}

*τ*/

*D*, where

*D*is the damping coefficient that captures the inertia of the system induced by the gating kinetics, as shown by FitzHugh [21], dimensional time is represented by

*τ*, and

*R*

_{1}is the passive resistance of the nonlinear current device of the circuit, represented as a tunnel diode by Nagumo

*et al.*[22]. Dimensionless coupling is defined as

*κ*=

*D*/

*R*

_{1}

*K*. The dimensionless equations governing the dynamics of the cell in the (

*i*,

*j*)th position are as follows: 2.8aand 2.8bHere

*I*is the input current applied to cell (

_{i,j}*i, j*). To render the remaining parameters dimensionless, define where

*R*

_{1}is the passive resistance of the nonlinear element,

*C*is the capacitance of cell (

_{i,j}*i*,

*j*) and

*D*is the damping coefficient. The quantities

*R*and

*D*are incorporated into the excitable element depicted in figure 2

*d*. The other parameters were scaled in accordance with the derivation given by Keener & Sneyd [28], as detailed in appendix A.

#### 2.2.4. Initial conditions and activation clusters

The initial conditions for each cell correspond to the resting membrane potential and are given by the real solution to the following simultaneous equations:
2.9and
2.10A perturbation *Δ**v*, representing a short-lasting influx of charge, is applied to the cell at the centre of the lattice, displacing it from its stationary point. The cells that become excited as a consequence of this initial perturbation constitute an activation cluster. The ratio of the number of cells in the cluster to the total number of cells in the simulated lattice is used as a measure of the strength of activation.

#### 2.2.5. Spatial structuring

Both homogeneously coupled and heterogeneously coupled cell networks are considered. In the spatially homogeneous case, all couplings are equal, i.e. *κ*_{i}_{,j} = *κ* ∀(*i*, *j*), whereas in the spatially heterogeneous case, the cell–cell couplings *κ*_{ij} are allowed to vary with *i*,*j*. Two further subcases can be distinguished: (i) symmetric coupling, i.e. *κ*_{i,j} = *κ*_{j,i} ∀(*i*, *j*) and (ii) asymmetric coupling, where *κ*_{i,j} ≠ *κ*_{j,i}.

##### 2.2.5.1. Symmetrical coupling: the Bernoulli lattice

To simulate symmetrical coupling, the value *κ*_{i,j} = *κ*_{j,i} is determined probabilistically by
2.11where the limiting case *p* = 1 represents spatial homogeneity. We refer to the lattice defined by equation (2.11) as the Bernoulli lattice. Heterogeneity in the Bernoulli lattice is modulated by varying the amount of connectivity in the lattice, expressed by the parameter *p* representing the probability of a connection existing between any two given cells.

##### 2.2.5.2. Symmetrical coupling: uniformly distributed

A second type of symmetrical heterogeneity is introduced using a lattice with connections drawn from uniform distributions. In particular, the intercellular couplings are modelled as independent uniform variates on [0.5,10].

##### 2.2.5.3. Spatial correlation in the coupling structure

Distance to the stimulated cell is expressed using the *ℓ*_{1} (‘city-block’) measure; electronic supplementary material, figure S1, labels cells according to their *ℓ*_{1}-distance to the central cell on a 5 × 5 lattice. The correlation function *C _{r}*, where

*r*is the

*ℓ*

_{1}-distance, is defined as follows: 2.12where denotes the grand mean of coupling over the network and is the mean coupling of cell

*i*.

##### 2.2.5.4. Asymmetrical coupling: cell capacitances

Asymmetrical coupling is introduced by allowing variation in cell capacitances (with symmetrical gap junctional conductances). A larger cell has a larger capacitance, and as a result will be more weakly coupled to its neighbours since *K _{i}*

_{,j}= (

*C*

_{i}R_{i}_{,j})

^{–1}will be smaller (figure 2

*e*). The model assumes that the cell capacitance values are independently lognormally distributed; calculation of realistic parameter values is detailed in §2.1.2. A resistance value is chosen from within a window of values that permit global excitability.

The value of the scaled parameter *ε* is proportional to the capacitance of the cell. Accordingly, in simulations using day 18 data, the value of *ε* is normalized by the modal point of the 18-day distribution. Similarly, *ε* is normalized by the modal point of the 15-day distribution in day 15 simulations.

##### 2.2.5.5. Variability in resting membrane potential

A normal distribution was fit to the normalized resting membrane potential data. The 18-day data serve as a reference point. The non-dimensional resting membrane potential for 18-day pregnant data is normally distributed with mean 1 and s.d. 0.1615, whereas the non-dimensional resting membrane potential for 15-day pregnant data is normally distributed with mean 1.046 and s.d. 0.08737.

Equations (2.9) and (2.10) can be solved simultaneously to find an expression for the value of *v* in the steady state. Retaining only linear terms, we find that the scaled resting membrane potential is proportional to *v*_{0} – *γw*_{0}.

A value *v _{i}*

_{,j}is sampled from the resting membrane potential distribution for each cell (

*i*,

*j*) in the network. The dependent parameters

*γ*and

*w*

_{0}are then multiplied by

*v*

_{i}_{,j}to introduce variability into the model. Introducing variability through the parameter

*v*

_{0}produces non-sensical results.

Multiplying *γ* and *w*_{0} by *v _{i}*

_{,j}alters the phase portrait for each cell. We therefore examine the effect of the variability in resting membrane potential (and so gestation period) on the minimum perturbation required to excite an isolated cell. Simulations were run with the only asymmetry being the variation in resting membrane potential, and also with an asymmetry in the cell capacitances. Relative cluster sizes were recorded for a range of coupling strengths and initial perturbations with varying connectivity

*p*.

An increase in variability in resting membrane potential is seen as pregnancy progresses. Therefore, by assuming a linear relationship between day 15 and day 18 data, we extrapolate resting membrane potential variability for time points outside this range to examine the trend of excitability throughout pregnancy.

##### 2.2.5.6. Pacemaker cells

Activity waves can be initiated by external stimuli (for example, oxytocin) or by intrinsic activity of spontaneously active cells. Such cells can be obtained by adjusting the parameters to shift the position of its null isoclines. The pacemaker cell parameter values are given in table 2. The straight line null isocline passes directly through the local minimum of the cubic null isocline to create an unstable fixed point at the rest state of the cell, as indicated in figure 1*d*. In effect, the cell continually re-excites itself. The cells are connected symmetrically with a constant coupling value.

The models used in this study are summarized in table 3.

## 3. Results

To investigate how the spatial heterogeneity affects the global excitability of the network, the spatially homogeneous case is studied first, followed by the heterogeneous coupling case.

### 3.1. Fully connected square lattice

Simulations of a spatially homogeneous square lattice indicate that for each perturbation there is a window of coupling strengths in which global excitation of the lattice is possible. For instance, an initial perturbation of 1 in a fully connected lattice of size 25 × 25, produces global excitation for coupling strengths (*κ*) between 0.76 and 5.12. Below 0.76, the excitation does not spread beyond the cell to which the stimulus was applied. Above 5.12, no cell is able to become excited, including the perturbed cell. The surrounding cells act as a current sink; when the coupling is too strong, this effect prevents excitation. The symbols *κ*_{min} and *κ*_{max} denote the boundaries of the window of excitation; for the standard perturbation *Δ**v* = 1 these values are 0.76 and 5.12.

The ranges of perturbation and coupling values that permitted global excitation are shown in figure 3 for a 25 × 25 lattice. Simulations run for smaller and larger lattices produce qualitatively similar excitation curves.

### 3.2. Symmetric coupling in the Bernoulli lattice

Simulations were run with a lattice of 25 × 25 cells and the coupling strengths were chosen in the range from 1 to 2.5 (i.e. within the window of global excitation for an initial perturbation of 1). Initial perturbations took the following values: 1, 1.5, 2 and 2.5. As the connectivity of the lattice (*p*) increases, the relative cluster size increases in a sigmoidal fashion, as shown in figure 4. The response curve does not substantially vary with lattice size, coupling value or initial perturbation (see figure 4 and electronic supplementary material, figure S2), which suggests that connectivity is the dominant factor governing the excitability of the system.

Coupling values were chosen from outside the window of global excitation, and the effect of spatial heterogeneity in coupling was investigated by varying lattice connectivity as before. The behaviour for coupling values greater than or equal to the maximum coupling value *κ*_{max} is shown in figure 5*b*–*d* for initial perturbation *Δ**v* = 1. The behaviour for the minimum coupling value *κ*_{min} is shown in figure 5*a*. In both the cases, the sigmoidal curve was replaced by a negatively skewed bell curve, indicating that even when global excitation is precluded in the fully connected lattice, it can be attained in a partially connected system. Connectivity attains an optimum in the range *p* = 0.6–0.8. At these probabilities, either more than 90 per cent of the cells or fewer than 10 per cent of cells become excited. The averaging of the cluster sizes over 100 replicates underlies the decline in cluster size seen in the bell-shaped curve. At *p* = 1, this all-or-nothing behaviour is absent. Lattice size had no effect on this result (see the electronic supplementary material, figure S3).

### 3.3. Global transitions

Simulations around boundary points *κ*_{min} and *κ*_{max} were performed to elucidate the transition from quiescence to global excitation, at varying levels of connectivity, *viz. p* = 1, *p* = 0.7, *p* = 0.5 and *p* = 0.3. Initially, simulations were run at *κ*_{max} with perturbations from 0.5 to 1.5 at probabilities of 1, 0.7, 0.5 and 0.3 (figure 6*a*). In the fully connected lattice (*p* = 1), no excitation occurs for *Δ**v* < 1, whereas global excitation is attained for *Δ**v* > 1. At reduced connectivity, the lattice exhibits global excitation even at perturbations *Δ**v* < 1. The reduced connectivity allows a more gradual increase in cluster size with increased stimulation. This effect was especially pronounced at *p* = 0.7, an almost exactly optimal level of connectivity.

A similar phenomenon occurs near the lower bound (figure 6*b*). The fully connected lattice is unexcitable when the coupling is below *κ*_{min}, whereas global excitation can occur at couplings above this value. At reduced levels of connectivity, coupling values below the lower bound are associated with activation clusters of increasing size. The strength of the perturbation, however, has little effect on the cluster size (see the electronic supplementary material, figure S4), indicating that connectivity *per se* is the dominant factor governing the transition to global excitability.

### 3.4. Uniformly distributed coupling strengths

The sum and standard deviation of the coupling strengths of the four cells surrounding the excited cell were examined for the fully connected lattice with couplings sampled from a uniform distribution between 0.5 and 10. For a combined coupling value over the four neighbouring cells (*κ*_{T}) less than 21.0, the cluster always exhibits global excitation, with every cell becoming excited. For *κ*_{T} greater than 22.2, the couplings around the perturbed cell do not allow any neighbouring cell to become excited. This points to a threshold for global excitation in a fully connected system with coupling sizes sampled from a uniform distribution. For *κ*_{T} between these two values, an increase in the standard deviation between the coupling strengths is correlated with a reduced global threshold for excitation (see the electronic supplementary material, table S2). The ability to achieve global excitation does not depend on the combined strength of the couplings to the neighbouring cells. The relationship between standard deviation and global excitation around the upper bound of coupling strengths (§3.1) for *Δ**v* = 1 is shown in table 4. The standard deviation has a more pronounced effect on the cluster size than *κ*_{T}. At constant *κ*_{T}, there is a standard deviation threshold beyond which the system cannot achieve global excitation.

### 3.5. Uniformly distributed coupling strengths: Bernoulli lattice

The two sources of heterogeneity were combined to give a lattice with coupling values drawn from a mixed distribution with a finite probability at the value zero and a uniform distribution over non-zero values. The isotropic lattice was studied for various values of the initial perturbation (figure 7*a*). The negatively skewed bell-shaped curve was observed with a fixed coupling value for the smallest perturbation. Compared with the Bernoulli lattice, larger clusters are observed, at higher connectivity values. Again, the declining slope of the bell-shaped curve for a perturbation of 1 is associated with all-or-nothing behaviour. The transition from quiescence to global excitability as a function of *p*, at fixed *Δ**v*, is gradual, in contrast to the Bernoulli lattice (figure 7*b*).

#### 3.5.1. Role of spatial correlation in the coupling structure

The spatial correlation function exhibits a striking difference between networks in which more than 90 per cent of the cells become excited, and those in which fewer than 10 per cent of the cells become excited. The latter are characterized by a strong correlation between the stimulated cell and the *r* = 1 neighbouring cells. The strong local coupling strengths act as a sink, preventing the current from dissipating throughout the network (as illustrated in electronic supplementary material, figure S5*a* where couplings are strong between the central cell and its *r* = 1 neighbours). In contrast, a low degree of spatial correlation between the stimulated cell and its *r* = 1 neighbour is associated with networks in which 90 per cent of cells become excited. An example of such a network is given in the electronic supplementary material, figure S5*b*. As the heterogeneity decreases from *p* = 0.6 to *p* = 0.9 (figure 8*a*–*d*) this effect becomes less pronounced. At a connectivity of *p* = 1, this difference in correlation is not observed (figure 8*e*).

### 3.6. Asymmetrical coupling: cell capacitance

The couplings between cells can be adjusted to represent differences in cell size as explained in §2.2.5. Realistic distributions were derived from data obtained in 15-day pregnant and 18-day pregnant mice, shown in electronic supplementary material, figure S6. Figure 9 shows a model with scaled *ε* values, a constant gap junctional resistance and lognormally distributed capacitance values. At perturbations of 1.5, 2.0 and 2.5, a sigmoidal curve is observed, with an increased probability of cell-to-cell coupling corresponding to a larger cluster size. For *Δ**v* = 1, a bell-shaped curve is observed, with slightly reduced connectivity resulting in a larger relative cluster size. The graph shown is for the day-18 simulations. There was no discernible difference between day-15 and day-18 simulations. It appears that this variation between gestation days is not sufficient to remove the cells from the excitable range displayed in figure 3. In contrast to the symmetrically coupled case discussed earlier (figure 4), the curve is not sigmoidal for an initial perturbation of 1, suggesting that an asymmetrically coupled system requires a larger initial perturbation before full excitation can be achieved at full connectivity.

### 3.7. Variation in resting membrane potential

If the resting membrane potential is allowed to vary between cells, the pattern of spatial heterogeneity enhancing excitability persists. Figure 10*a* indicates the value of *Δ**v* required to excite the cell for a given resting membrane potential mean and standard deviation. The means and standard deviations for other gestation days were determined by interpolation and extrapolation, assuming a linear relationship between the 15-day and 18-day data. With the change in resting membrane potential, the cells' phase portraits change over time, resulting in a smaller *Δ**v*-value required for excitation of an individual cell as the system approaches parturition.

Figure 10*b* illustrates the relationship between relative cluster size and connectivity for systems with pre-day 15, day 15, day 18 and post-day 18 resting membrane potential variation. We see that a system with resting membrane potential variation only, at a constant coupling strength (*κ* = 1), fails to achieve full excitation. However, increased gestational age correlates with increased excitability.

Figure 11 displays the results of variation in resting membrane potential and cell capacitance, for all four gestation time point combinations. Capacitance variation endows some cells with a greater ability to excite the network, since larger cells can accommodate the current sink effect exerted by its neighbours and act as a buffer, enabling all cells to become excited. It appears that excitability is optimal with larger perturbations in the case combining capacitance and resting membrane potential variability sampled from the 18-day distribution. In contrast to the symmetrically coupled case (figure 4), an asymmetrically coupled system requires a larger initial perturbation before full excitation can be achieved at full excitability.

### 3.8. Pacemaker cells

In networks with a pacemaker cell, the frequency of oscillations of the pacemaker is affected by connectivity. The current sink represented by neighbouring cells can reduce the frequency of the pacemaker cell's free-running oscillation or even cause the cell to cease cycling altogether. The probability of such a complete cessation of cycling increases with connectivity (figure 12*a*). With the exception of very small coupling values *κ*, at full connectivity the pacemaker cell is not able to maintain its excitability and subsequently none of the surrounding cells are able to become excited. At small values of *κ*, the pacemaker maintains its excitability and pacemaker frequency is hardly affected by increased probability of connection. Heterogeneity of the spatial coupling structure thus appears to be an important modulator of excitability.

The relationship between connectivity and the frequency of oscillations of the pacemaker cell is shown in figure 12*b*, calculated as an average of simulations in which the pacemaker continues to cycle at a finite frequency. As the connectivity *p* increases from 0 to 1, the frequency of oscillations of the pacemaker cell decreases, owing to the drain of current to the surrounding connections. At high coupling strength values of *κ* = 5 and *κ* = 10, connectivity has no effect on frequency since the only way for a pacemaker to retain a finite frequency at these coupling strengths is to be essentially isolated in the network. This is illustrated in electronic supplementary material, figure S7 that shows the topology of a sample pacemaker network with finite frequency, and the topology of a sample network where the pacemaker's frequency decays.

The key results presented in this paper are summarized in table 3.

## 4. Discussion

It is well known that an increase in the number of gap junctions in myometrial tissue is indicative of the onset of labour [1]. While the number of connections increases, little is known about the strength and local structure of these connections. The main finding of this study is that heterogeneity of the coupling structure of a network of excitable elements allows the system to respond in a graded manner to a wide range of stimuli. For instance, reduced connectivity results in a gradual increase in cluster size in response to stimulation, resulting in a smooth transition between an unexcited and globally excited state, which enhances the scope for precise regulation. Physiological evidence suggests that the transition from quiescence to excitation is gradual [16]; this may be explained in part by the results presented here.

Relatively little is known about the effect of heterogeneity in myocyte connections on the excitability of the network as pregnancy progresses. Our analysis of the pregnant mouse data confirms the importance of heterogeneity and demonstrates that the network evolves towards global excitability as labour approaches.

A homogeneous lattice (fully connected with identical couplings) cannot achieve global excitation when the coupling is too strong: neighbouring cells can too readily absorb the membrane charge associated with the stimulus. When heterogeneity is introduced, whether by removing selected couplings at random or by varying the coupling values at random, the ability to achieve global excitation is restored. An intuitive explanation is that heterogeneity creates local ‘pockets’ of excitability where the initially stimulated cell is less hampered by the charge drain imposed by its neighbours. Moreover, heterogeneity smoothes out the sudden transition between an unexcited network and a globally excited one (§3.3). This is consistent with the finding that in the days and hours leading up to labour, contractions spread further throughout the uterus [29]. Finally, spatial heterogeneity regulates the ability of pacemaker cells to drive the network, again in a graded fashion.

Simulations on a fully connected square lattice (§3.1) indicate that there exist minimum and maximum coupling values between which global excitation of the network is possible. These thresholds were found to be dependent on the initial perturbation. Whereas the homogeneous lattice cannot exhibit global excitation, the heterogeneous lattice formed by removing couplings is able to exhibit excitation. In the heterogeneous system, the wave of excitation generated by the central cell has to overcome a lower threshold to excite neighbouring cells. Therefore, the system can achieve global excitation even with a smaller inciting stimulus.

Investigation of the variability of the couplings surrounding the excited cell in a lattice with couplings drawn from a uniform distribution indicates that global excitation is affected by the variance of the coupling strengths to neighbouring cells and is not solely dependent upon total coupling. The effect of coupling variance displayed a lower and upper threshold for global excitation of the network. Therefore, a variance within these bounds is optimal for network excitation. Outside this region, excitability does not spread across the whole network. It appears that when the variance of connections is too high, a fragmented network results, which does not allow for propagation of excitability. Overall, it appears that variation of coupling allows for a larger cluster size at high connectivity. This result was supported by similar simulations run for an anisotropic lattice (see the electronic supplementary material, table S2), where the horizontal and vertical coupling strengths were drawn from distinct uniform distributions.

The importance of heterogeneity was demonstrated by the correlation between cells at increasing distance from the central cell, and the central cell itself. At full connectivity, there was no discernible difference in correlation. For networks with a lower probability of connection, a much larger correlation was observed in networks with fewer than 10 per cent of cells becoming excited when compared with networks with more than 90 per cent of cells becoming excited. This effect was more pronounced in lower connectivities, reaffirming the importance of heterogeneity. This supports the notion of excitation pockets, where a cluster of cells that are too well-connected act as a current sink and prevent the excitation from spreading further through the network. The main determinant of global excitation is the local spatial correlation around the excited cell.

Data for resting membrane potential in mice at 15-day and 18-day gestational age indicate an increase in depolarization of the cell membrane with increased gestation length, confirming the findings of Lodge & Sproat [13]. The statistical variation in resting membrane potential among cells also increases with gestational age. With resting membrane potential as the only form of variability in the model, full excitation could not be reached. However, when the resting membrane potential's variability is combined with capacitance variation, the optimum excitability at intermediate connectivity is restored.

Results were consistent over a range of lattice sizes; from 9 × 9 to 25 × 25. Larger lattices were not investigated due to computational restrictions. Simulations were also run for hexagonal lattices to investigate effects of network structure on excitability. Similar qualitative results were obtained. However, the hexagonal lattice displayed a wider window of global excitability when compared with a square lattice with the same number of nodes. An increase in spatial heterogeneity enables the system to achieve global excitation more easily. Networks with stochastically varying node degrees will be investigated in future research.

The resting membrane potential values were taken from recordings of isolated myometrial cells. However, a system of connected cells settles into a different resting states. To verify that the experimental data were representative of networked cells, simulations were run in which the cells where allowed to equilibrate, and differences were found to be negligible.

Preliminary observational data (shown in the electronic supplementary material, section 10.1 and figure S8) suggest that an increase in heterogeneity increases the frequency of the oscillations of the pacemaker. This was confirmed by our simulations in which the pacemaker cell's ability to drive the system is modulated by connectivity. With increased heterogeneity, the drain of current to surrounding cells has a smaller effect, causing the frequency of oscillations to increase. In addition, a well-connected system makes it more probable that the pacemaker is not able to remain active, and its oscillations stop altogether. Maximal excitability is associated with a moderate degree of spatial heterogeneity. This would suggest that a way in which a cell can develop pacemaker activity is to downregulate its gap junctions to partially isolate itself in the network. The suggestion of cell isolation is supported by preliminary experimental data (not shown) in which individual cells can oscillate but still participate in a global action potential.

Propagation waves in excitable media have been examined extensively using cellular automata [30–34]. This approach was then extended to show how the cellular automata models could be matched to systems such as the FitzHugh–Nagumo set of equations [35,36]. A cellular automata model also exists that describes the uterine network using Hodgkin–Huxley physiology [37]. Here, we have opted to use the FitzHugh–Nagumo model to represent the excitation and recovery of a myometrial cell. Since the primary aim of the present study was to explore the consequences of spatial heterogeneity *per se*, we have used a two-dimensional minimal model to represent the excitable element without being too computationally expensive. This FitzHugh–Nagumo model captures the qualitative dynamics of propagation through excitable systems. It does, however, still have its limitations. The model considers only one activation variable and one recovery variable. The minimal model can be replaced by models that take into account individual currents and the intracellular calcium stores [38–40]. Future work will focus on these more detailed models to more accurately represent the calcium influx resulting in a contraction. In addition, a three-dimensional system will be considered to model the uterine network more closely.

Miyoshi *et al.* [41] demonstrated that gap junctional conductance between coupled cells is dependent on the trans-junctional voltage, which adds a nonlinear effect that may play an important role in the transition from the globally quiescent state to the globally excitable state. This effect has not been taken into account in this study and will be the subject of future research.

In summary, the mathematical model used here indicates that spatial heterogeneity may serve as an important modulator of excitability in uterine muscle. Spatial heterogeneity of cell-to-cell connections promotes an increase in the excitability of the network, and the ability of a network to become fully excited is governed predominantly by the local connection structure. In addition, heterogeneity in both cell capacitance and resting membrane potential also plays a role. Similarly, heterogeneity allows a pacemaker cell to drive the system.

Shifts in this heterogeneity may be a significant factor in the regulation of myometrial excitability as pregnancy progresses from conception to parturition. Pre-term labour may be associated with a premature development of spatial heterogeneity. Mapping of spatial heterogeneity may prove to be a diagnostic tool to monitor the development of excitability throughout pregnancy.

## Acknowledgements

R.E.S. gratefully acknowledges funding for a studentship from the Engineering and Physical Sciences Research Council through the MOAC Doctoral Training Centre at the University of Warwick.

## Appendix A

We follow the derivation provided by Keener & Sneyd [28]. Consider the circuit given in figure 2*d* showing a simplified model of the cell membrane. From Kirchoff's laws, one obtains
A1and
A2where *I*_{0} is the applied external current, *i* is the current through the resistor–inductor, *V* = *V _{i}* –

*V*is the membrane potential,

_{e}*R*is the resistance and

*V*

_{0}is the potential gain across the battery, which forms part of the excitable element representing the recovery current in figure 2

*d*.

*D*is the damping coefficient, which captures the inertia of the system induced by the gating kinetics, as shown by FitzHugh in 1961 [21].

*R*and

*D*are incorporated into the ‘excitable element’. Here,

*τ*represents dimensional time. The function

*F*(

*V*) is a cubic with three zeros: the smallest

*V*= 0 and largest

*V*=

*V*

_{1}are stable solutions of d

*V*/d

*τ*= –

*F*(

*V*). The passive resistance of the nonlinear element (defined as a tunnel diode by Nagumo [22]) is

*R*

_{1}= 1/

*F*′(0).

The equations are rendered dimensionless as follows. Define *v* = *V*/*V*_{1}, *w* = *R*_{1}*i*/*V*_{1}, *f*(*v*) = –*R*_{1}*F*(*V*_{1}*v*)/*V*_{1} and *t* = *R*_{1}*τ*/*L*. Equations (A 1, A 2) can then be rewritten as follows:
A3and
A4where , *w*_{0} = *R*_{1}*I*_{0}/*V*_{1}, *v*_{0} = *V*_{0}/*V*_{1} and *γ* = *R*/*R*_{1}. The function *f*(*v*) is a cubic and can be written as follows:
A5

- Received May 20, 2013.
- Accepted June 19, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.