## Abstract

The smooth muscle cells of the uterus contract in unison during delivery. These cells achieve coordinated activity via electrical connections called gap junctions which consist of aggregated connexin proteins such as connexin43 and connexin45. The density of gap junctions governs the excitability of the myometrium (among other factors). An increase in gap junction density occurs immediately prior to parturition. We extend a mathematical model of the myometrium by incorporating the voltage-dependence of gap junctions that has been demonstrated in the experimental literature. Two functional subtypes exist, corresponding to systems with predominantly connexin43 and predominantly connexin45, respectively. Our simulation results indicate that the gap junction protein connexin45 acts as a negative modulator of uterine excitability, and hence, activity. A network with a higher proportion of connexin45 relative to connexin43 is unable to excite every cell. Connexin45 has much more rapid gating kinetics than connexin43 which we show limits the maximum duration of a local burst of activity. We propose that this effect regulates the degree of synchronous excitation attained during a contraction. Our results support the hypothesis that as labour approaches, connexin45 is downregulated to allow action potentials to spread more readily through the myometrium.

## 1. Introduction

The myometrium, which makes up the bulk of the uterine wall in all species [1] (electronic supplementary material, figures S1*a* and S1*b*), is composed of interconnected smooth muscle cells, which form an excitable medium [2], i.e. a medium that is capable of propagating signals over relatively long distances with virtually no damping. In the myometrium, these propagating signals take the form of action potentials which trigger phasic contractions during labour [3]. These involve the rhythmic contraction and relaxation of the smooth muscle cells, giving rise to synchronized travelling waves throughout the uterus [4]. In humans, the phasic contractions typically last for around a minute [5], whereas in rats, they last for approximately 68 s at term [6].

Throughout most of pregnancy, the myometrium remains quiescent, allowing the fetus to develop undisturbed [7]. The myometrium undergoes major changes in the days leading up to labour [7]. In particular, as parturition approaches, the uterus becomes more excitable through an activation phase which involves molecular changes that lead to a decrease in time between subsequent contractions [8]. In humans, this interval decreases from approximately 30 min to approximately 2 min as labour progresses [5,8], whereas in rats, contractions occur at 2.3 min intervals at gestational day 19 and at 46 s intervals in active labour [6].

A key condition for the transition to excitability at the whole-organ level is a change in connectivity between the myometrial smooth muscle cells [2]. The smooth muscle cells in the myometrium are connected through a series of intercellular channels that allow the passage of ions and small molecules from one cell to another, each channel being made up of two components, one contributed by each cell [9].

As pregnancy progresses, the density of gap junctions on the cell-to-cell contact surfaces increases, which is thought to be one of many factors that instigate labour [2]. In the non-pregnant uterus of both rats and humans, gap junctions are either absent or present at low density [2,10,11]. At the time of delivery, gap junctions are present at a density of approximately 1000 per cell in human tissue [12].

The ability of a myometrial network to achieve global excitation is strongly dependent on spatial heterogeneity of cell-to-cell connectivity [13]. Using a mathematical model of the myometrium based on FitzHugh–Nagumo dynamics [14,15], we previously showed that a moderate degree of spatial heterogeneity allows excitation to spread throughout a network even if the same network would be quiescent when fully and uniformly connected. However, the coupling between cells was represented as a passive resistor, whereas, in reality, voltage-dependent gating is a well-established property of vertebrate gap junction channels [16–19]. This transjunctional voltage-dependence is directly related to the connexin composition of the channel. Individual connexins exhibit different responses as a function of the potential difference between neighbouring cells; for example, connexin43 (Cx43) is relatively insensitive to the voltage difference, displaying only modest changes in conductance as the transjunctional voltage is altered [20]. By contrast, connexin45 (Cx45) and connexin40 are highly sensitive to transjunctional voltage differences [20–22]. The differences in connexin protein response to transjunctional voltage are thought to be underpinned by different localizations of charged residues along the inner lining of the ion-conducting pore as well as different phosphorylation states [23,24].

Miyoshi *et al.* [25] recorded gap-junctional currents between isolated pairs of rat myometrial cells, using the double-whole-cell voltage-clamp configuration (electronic supplementary material, figure S3). The cell pairs were found to exhibit one of two distinct relationships between conductance and voltage-dependence which occurred in equal proportions (figure 1*a*). The so-called type I channels displayed a gradual increase from a normalized conductance of 0.2–1 with decreasing voltage difference. At small transjunctional voltages, there was very little change in conductance. By contrast, the ‘type II’ channels showed a much sharper increase in conductance with decreasing voltage difference. A notable voltage-dependent change in conductance was observed even at small transjunctional voltages. The authors observed that the type I voltage-dependence relationship fitted with the pattern expected for a tissue with predominantly Cx43 junction proteins. In addition, the type II junction voltage-dependence appeared to be notably distinct from the connexin patterns observed before. The dependence of conductance on the transjunctional voltage difference for type II is more consistent with a Cx45 junction protein, a hypothesis that was further supported by the findings by Moreno *et al.* [26].

Here, we investigate how the spread of excitation through a network is affected by incorporating voltage-dependent cell-to-cell couplings in the mathematical model. We adapt the voltage-dependent conductance relationships for both gap junction types identified by Miyoshi *et al.* [25]. The main qualitative effect is that conductance increases as the potential difference between the cells decreases. To isolate this ‘cut-off’ effect, we initially model the voltage-dependence relationship as a simple step function. We subsequently compare our findings with more realistic voltage-conductance relationships in order to verify that the dominant qualitative effect is the ‘cut-off’ point of activity.

Finally, we incorporate gating kinetics. We propose that the type II gap junctions are not able to excite the network by themselves; a substantial proportion of type I gap junctions is needed to spread activity effectively. The gating kinetics is much more rapid for type II gap junctions in comparison with type I junctions. The faster kinetics prevents the longer depolarization propagating throughout the network, and so limits the excitability of the network. Accordingly, we suggest that relative proportions of type I and type II gap junctions are crucial for the modulation of activity in the excitable myometrial network.

Analysis of gene expression in rat models suggests that the expression of Cx45 is downregulated in labour. We do not observe a change in the expression of Cx43 in the transition from pregnancy to labour. We therefore hypothesize that the relative proportion of Cx43 and Cx45 is crucial; Cx45 acts as a brake on activity in the uterus, and its downregulation may be an important step in the initiation of labour.

## 2. Methods

### 2.1. Mathematical modelling

We further develop a mathematical model of the myometrial smooth muscle network [13]. The network is constructed as a square lattice in which each node is connected to four others. Each simulated cell obeys FitzHugh–Nagumo dynamics. The FitzHugh–Nagumo model is a two-variable system that describes the excitation and relaxation of an excitable cell [27]. The model is an example of an excitation–relaxation oscillator owing to its threshold properties in which a sufficiently large input of current is required for the cell to exhibit an excursion corresponding to excitation and relaxation. However, parameters can be chosen such that cells are autorhythmic [13], that is, the cells continually oscillate through the excitation–relaxation cycle without any external stimulus. FitzHugh–Nagumo is a simplification of the Hodgkin–Huxley model for spiking neurons [28,29] and, as such, is often used in modelling excitable systems such as contractions of cardiac myocytes [30–32].

Each node in the network is considered to be an excitation–relaxation oscillator with a fast (‘excitation’) state variable *v* that corresponds to the transmembrane potential and a slow (‘recovery’) variable *w* that corresponds to gating kinetics which repolarize the excited cell. The following dimensionless ordinary differential equations describe this two-variable model
2.1aand
2.1bwhere *I* is the input current and *B*, *α*, *γ* , *w*_{0}, *v*_{0}, and *ε* are positive parameters. The values shown in table 1 are used in all simulations to generate a network of cells which require a suprathreshold current input to undergo an oscillation.

The equations were solved using a Livermore solver for ordinary differential equations approach through the NDSolve function in Mathematica. Results were numerically stable under variation in step size.

Each simulated cell is coupled to its four neighbours through resistors representing gap junctions (electronic supplementary material, figure S2*a*). A coupling constant *K* is assigned to each gap junction and defined as *K _{ij}* = (

*C*)

_{i}R_{ij}^{−1}, where

*i*and

*j*are two neighbouring cells,

*C*is the capacitance of cell

_{i}*i*, and

*R*is the resistance between the two cells. The coupling

_{ij}*K*is therefore a rate constant, and is rendered dimensionless as follows. The FitzHugh–Nagumo model defines dimensionless time,

*t*, as

*t*=

*R*

_{1}/

*τD*, where

*D*is a damping coefficient that captures the inertia of the system owing to the gating kinetics, as shown by FitzHugh [14],

*τ*is dimensional time and

*R*

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

*et al.*[15]. Accordingly, dimensionless coupling is defined as

*D*/

*R*

_{1}

*K*and denoted by

*κ*. The scaling of the remaining parameters is detailed in the appendix and follows the derivation given by Keener & Sneyd [27].

Miyoshi *et al.* [25] identified two distinct relationships between the conductance of smooth muscle gap junctions and the voltage difference between two neighbouring myocytes (figure 1*a*). We adapt these relationships to represent voltage-dependent gap junctions in terms of the coupling strength *κ*. Spatial heterogeneity is added to the system by introducing a probability *p* of a connection existing between any two adjacent cells, as previously described [13]. This parameter is referred to as the ‘connectivity’ of the network. The relative cluster size serves as a measure of the proportion of cells in the network that become excited: if all cells are excited, the relative cluster size equals 1, whereas, if no cells become excited, it is 0.

#### 2.1.1. Step function model

The simplest way to represent the voltage-dependence is by means of a step function (figure 1*b*). For cells with a small voltage difference, the normalized gap-junctional conductance is set to be 1, whereas for cells with a large voltage difference, the normalized gap-junctional conductance is 0.2 (the base value determined by Miyoshi *et al.* [25]). We use the term *gap-junctional conductance bandwidth* to mean the range of voltage differences for which the gap-junctional conductance equals 1. The voltage differences are scaled from unit-bearing values to the dimensionless FitzHugh–Nagumo model values. To facilitate numerical analysis, the step function is smoothed, as follows (electronic supplementary material, figure S4*a,b*)
2.2where the indices *i* and *j* represent adjacent cells, *λ* is the scaling required to convert between physiological values and the FitzHugh–Nagumo model values, and *ψ* is the gap-junctional conductance bandwidth. The relative conductance is multiplied by the coupling strength *κ*. In the majority of simulations, *κ* = 1; however, we also present some results with *κ* = 0.76 to facilitate comparison with previous results in which we demonstrated the importance of heterogeneity [13].

Scaling is based on a typical rat myometrial cell in accordance with Miyoshi *et al.* [25]. The membrane potential of a rat myometrial cell increases from −45 to +10 mV during the upstroke of a full action potential, corresponding to an excursion of 55 mV. In the non-dimensionalized model, the size of the excursion is 3.52278. Accordingly, *λ* is taken to be 55/3.52278. To vary the width of the gap-junctional conductance band *ψ* is selected from the range 40 to 200 mV (twice the difference to zero, which lies at the midpoint). The value of *ψ* then dictates the voltage difference values over which the conductance of the gap junction equals 1. For example, *ψ* = 120 mV sets the conductance to 1 over a voltage difference of ±60 mV, as illustrated in figure 1*b*.

To check the model against our previous findings, the bandwidth is set to a sufficiently high value such that the relative conductance remains 1 over the entire voltage range of the action potential. The proportion of cells that become excited as a function of connectivity *p* reproduce our previous results; see the electronic supplementary material, figure S5 for more details. With a coupling strength *κ* = 0.76, we replicate our previous scenarios in which spatial heterogeneity proved to be crucial in exciting the network (electronic supplementary material, figure S6).

#### 2.1.2. Symmetrical Miyoshi model

Miyoshi *et al.* [25] proposed the following empirical equation for the conductance–voltage relationship
2.3where *g* is the normalized conductance of the gap junction, *V*_{j} is the voltage difference across the gap junction, *G*_{min} is the minimum value of the conductance, *V*_{h} is the half inactivation voltage (mV) at which the conductance is reduced by (1 − *G*_{min})/2, and *A* is a slope factor ((mV)^{–1}). The authors identified two distinct types of voltage-dependent gap junctions with the parameters given in table 2. The relationship between the voltage gradient and gap-junctional conductance is shown for both types in figure 1*a*.

To investigate the effect of the width of the gap-junctional conductance band, we model the type I relationship by means of the bandwidth parameter *ψ* that allows us to scale the type I relationship to resemble the type II relationship. We employ a symmetrical version by taking the absolute value of the voltage and using the positive values that Miyoshi *et al.* [25] obtained for the Boltzmann fit, as shown in the equation (2.4)
2.4where *λ* is the scaling required to convert between physiological values and FitzHugh–Nagumo values, and *ψ* is the required width of the gap-junctional conductance band. For instance, *ψ* = 110 mV will give type I, and *ψ* = 40 mV will give type II. The scaling parameter *λ* = 55/3.52278, as before. A pictorial representation of the symmetrical Miyoshi model [25] is shown in figure 1*c*.

We checked the model against our previous results by setting the bandwidth to a high value in order to attain a relative conductance of 1 over the action potential voltage range (electronic supplementary material, figure S7).

#### 2.1.3. Symmetrical Miyoshi model with gating kinetics

The dependence of conductance on trans-gap-junctional voltage as described in §2.1.2 is non-dynamic, i.e. the gating kinetics is sufficiently rapid to justify a quasi-instantaneous model. We now turn to more realistic gating kinetics. In the first instance, voltage-dependence of the time constant is ignored, giving rise to the following ordinary differential equation
2.5where *g* is the conductance of an individual gap junction, *τ*_{g} is the gap-junctional time constant, and is the steady-state conductance dependent on the transjunctional voltage *V*_{j} as considered in §§2.1.1 and 2.1.2. A complete excitation–relaxation cycle in the FitzHugh–Nagumo model lasts for 18.6849 dimensionless units (the time at which the cell is able to become re-excited) and we assume an action potential burst in the rat lasts for around 30 s. Therefore, to convert from dimensionless time to seconds, we multiply the dimensionless time by 30/18.6849.

We next consider a model in which the time constant *τ*_{g} is dependent on the potential difference between neighbouring cells. The voltage-dependence relationship for the time constants is derived from Miyoshi *et al.* [25] by considering the gap-junctional current response after 5 s step pulses of voltages in the range ±90 mV (figure 2*a*,*b*, respectively). The type I and type II gap junctions exhibit distinct relationships, given in figure 2. An exponential equation is fitted using the least-squares criterion to each current response in order to determine the time constant for each voltage difference. A Gaussian is then fitted using the least-squares criterion to the time constants giving voltage-dependent time constant relationships as follows:
2.6where *V*_{j} is the transjunctional voltage, *a* = 9.39726, *b* = 23.9757, and *c* = 0.60274 for type I, and *a* = 9.41999, *b* = 16.799, and *c* = 0.580013 for type II (figure 3). We are not able to obtain good estimates of the rate of decay near 0 mV, because the time constant here substantially exceeds the duration of the experiment.

### 2.2. Experimental methods

We employed rat models in order to investigate the changes in expression of Cx43 and Cx45 in the myometrium between gestational ages of 19.25 and 22 days. Detailed experimental methods are provided in the electronic supplementary material, §S1.1.

## 3. Results

The effect of varying the gap-junctional conductance bandwidth is studied first. We then examine cases in which the conductance relaxes to a steady state over a certain period of time. In more sophisticated simulations, the time constant is dependent on the voltage difference between neighbouring cells.

### 3.1. Minimal gap-junctional conductance bandwidth

Numerical simulations with the step function model indicate that there is a critical minimal gap-junctional conductance band below which global excitation cannot be achieved (figure 4*a*). For conductance bandwidths of 100, 80, 60, and 40 mV, full activity is observed with a fully connected lattice. However, once the conductance bandwidth drops lower than 40 mV, all activity disappears, and only the perturbed cell is capable of being excited. The transition from excitation to quiescence is quite steep and takes place over a very small range of bandwidth values: between 40 and 39 mV (figure 4*b*). Simulations with several lattice sizes indicate that this transition has the character of a discontinuous phase transition in the limit of infinite lattice size (see electronic supplementary material, figures S8 and S9 for further details).

The critical value of the bandwidth is located between the apparent bandwidth values of type I and type II gap junctions as described by Miyoshi *et al.* [25] (figure 5*a*). This strongly suggests that physiological significance can be attached to the qualitative difference between bandwidths for which no global excitation is achieved, and those bandwidths that permit global excitation.

To assess the hypothesis that the threshold bandwidth value is significant, we ran simulations in which cell-to-cell couplings were randomly assigned the bandwidth corresponding to type I or type II, using probabilities of intercell connections of 0.4, 0.6, 0.8, and 1. The relative cluster size increases with the proportion of type I gap junctions (figure 6), with a steep increase between ratios of type II to type I gap junctions of 0.4 and 0.6. This qualitative result is also found with lower connection probabilities. As the connectivity of the lattice is decreased, the transition values shift in favour of a larger proportion of type I gap junctions.

The qualitative results obtained with the step function model were repeated in a model which more closely represented the shape of the type I and type II [25] fits (see §2.1.2). There is a minimum gap-junctional conductance bandwidth beyond which no excitation can take place, and a clear transition between excited and quiescent states. The results of this simulation are detailed in electronic supplementary material, figures S10 and S11.

Network activity depends on the connectivity *p* and the ratio of type II to type I gap junctions. Figure 7 suggests that there is a defined set of conditions that promote network excitation. The majority of parameter combinations result in no activity. However, a network with high connectivity and a higher proportion of type I gap junctions relative to type II gap junctions maximizes the likelihood of near full activity.

### 3.2. The effect of gating kinetics

Here, we examine the symmetrical Miyoshi model with voltage-gated gap junctions using a fixed time constant. The excitability of a network with exclusively type I gap junctions is unaffected by the value of the time constant—there is always full activity at full connectivity (figure 8*a*). By contrast, the time constant is a key parameter in a network with type II gap junctions. With a small time constant, type II gap junctions are unable to propagate any activity, as before. However, as the time constant is increased, the relationship between connectivity and relative cluster size assumes one of two distinct forms (figure 8*b*). For time constants between 0.87 and 1.49 s (0.54 and 0.93 in dimensionless units), the relationship follows a bell-shaped curve, indicative of the importance of spatial heterogeneity in exciting the network. For time constants larger than 1.51 s (0.94 in dimensionless units), the relationship follows a sigmoidal curve in which full connectivity corresponds to full excitation. These results are summarized in table 3.

Comparisons are drawn between the time constant thresholds identified above, and the time constant relationship with voltage difference for type II gap junctions (figure 3*b*). For a voltage difference greater than ±50 mV, excitation in the network is not possible. For voltage differences less than ±35 mV, excitation is always possible. When the voltage difference is between ±35 and ±50 mV, spatial heterogeneity in connections is essential to obtain excitation in the network.

Finally, we examine a network with varying proportions of type I and type II gap junctions. As before, cell-to-cell couplings were randomly assigned the bandwidth corresponding to type I or type II. In addition, the voltage-dependent time constants corresponding to the chosen bandwidth are incorporated into the model using the Gaussian fit detailed in equation (2.6). Type II gap junctions remain unable to produce full network excitation even at full connectivity (figure 9). A small influence of spatial heterogeneity is apparent for networks with exclusively type II gap junctions. Around 20% of cells are able to become excited with connectivities in the range 0.6–0.8. Type I gap junctions are unaffected by the switch to voltage-dependent time constants. The cells in the network are still able to exhibit full activity with full connectivity. As the connectivity is reduced, so is the excitation.

### 3.3. Experimental results

Transcriptomics data suggest that expression of connexin45 (Cx45) in the inner myometrium is substantially lower on gestation day 22 (in labour) when compared with gestation day 19.25 (figure 10*a*). On the other hand, there does not appear to be a difference in connexin43 (Cx43) expression between these two points in time in either the inner or outer myometrium (figure 10*b*). Taken together, these data would suggest an increase of the Cx43 : Cx45 ratio, i.e. a relative enrichment of Cx43.

## 4. Discussion

Our findings suggest that the relative proportions of connexin43 (Cx43) and connexin45 (Cx45) are vital in initiating contractions in the myometrium during parturition. We build on the results obtained by Miyoshi *et al.* [25] to show that global excitation can be attained in a model network relying solely on Cx43 proteins, whereas this is not the case for a system using only Cx45. In addition, we develop a model in which the conductance of a gap junction is dynamic and has a voltage-dependent time constant, Cx45 having faster gating kinetics than Cx43. We demonstrated that this time constant is a key factor in determining whether global excitation can be achieved. At shorter time constants, heterogeneity in the density of gap junction connections plays a more prominent role in deciding whether activity can spread.

In rodents (as in most mammals, but not humans [33]), the myometrium can be divided into two distinct layers of smooth muscle: the inner circular layer and the outer longitudinal layer [34]. Our experimental data indicate downregulation of Cx45 expression in the rat as labour approaches, in contrast to Cx43 expression where no change at the end of pregnancy appears to be evident. Chan *et al.* [35] found similar results in human myometrium, analysing the RNA of myometrial samples from pregnant women (five from women in labour, and five from women at full term, but not in labour). These authors showed that Cx43 expression remains unchanged between the two sample groups, whereas Cx45 expression is significantly downregulated in the ‘in labour’ group in comparison with the ‘not-in-labour’ group (electronic supplementary material, figure S12). These results are in accordance with our data, which furthermore suggest that the significant changes found by Chan *et al.* [35] were owing to changes in expression taking place in the inner myometrium as opposed to the outer myometrium. Our model simulations show that an elevated proportion of Cx43 proteins is a key factor in modulating conductance. This accords well with results obtained by Döring *et al.* [36] and Tong *et al.* [37], who found that in a Cx43 knockout mouse, the lack of intracellular myometrial coupling results in weaker contractions and delayed labour. These various lines of evidence strongly suggest that the Cx43 : Cx45 ratio is instrumental in controlling the global excitability of the myometriocyte network.

The connexin protein Cx45 plays an important role in cardiac myocytes; it is found in the atrioventricular (AV) node and the adjoining His bundles [38,39]. While Cx45 is not essential for continued survival in the adult mouse [39], Cx45 knockout mice die from heart failure *in utero* [40,41]. Without Cx45, contractions initiated in the AV node are not coordinated with the contractions in the ventricles; Cx45 may play a similar role in the myometrium, i.e. to inhibit contractions until the development of the fetus has advanced sufficiently for delivery. The inviability of Cx45 knockout mouse embryos [40,41] precludes a test of this hypothesis *in vivo*. However, Cx45 knockout *in vitro* studies could be used to demonstrate the development of global network events.

Voltage-dependent gap junctions have previously been considered in mathematical models [42,43]. Baigent *et al.* [44,45] and Donnell *et al.* [46] considered a two-cell model coupled by a dynamic gap junction which resides in one of three states, giving insights into cell connections on a local scale, emphasizing channel kinetics and ionic flows in great detail. By contrast, the present model focuses on overall excitability of the network as a function of spatial heterogeneity of the functional properties of the gap junctions, rather than on the propagation of activity wavefronts and the spatial patterning of such fronts; these aspects have been considered elsewhere [4,32,47,48]. Whereas we considered the global activation threshold in an electrotonically coupled smooth muscle syncytium, a general argument regarding the all-or-none character of synaptically coupled networks, such as the central nervous system, was advanced by Ashby *et al.* [49].

This work has not considered heteromeric gap-junction channels with both Cx45 and Cx43 proteins. We assumed that the relationships given by Miyoshi *et al.* [25] correspond to homomeric gap junction channels. Martinez *et al.* [50] demonstrated that these two gap junction proteins are coexpressed in cardiac cells, and that homomeric channels of either Cx45 or Cx43 have different unitary conductances. A heteromeric gap-junction channel might exhibit non-symmetrical behaviour; this will require further experimental studies in order to inform our model.

We do not report an increase in Cx43 expression in the rat between day 19.25 and day 22. This finding is in keeping with the literature in that the rise in Cx43 occurs prior to parturition (4 days before term in the rat). Cx43 proteins are stored in vesicles and only trafficked to the membrane gap junctions when labour commences, so causing an increase in functional expression in labour, but not an increase in gene expression in the tissue [11,51–55]. Studies using immunohistochemistry and Western blot analysis have demonstrated the functional increase in Cx43 expression at term as well as a decrease in Cx45 [56,57]. Albrecht *et al.* [56] suggested that the ratio of Cx43 to Cx45 might be crucial for increased network connectivity associated with the initiation of labour, which is in agreement with our finding that relative numbers of Cx43 and Cx45 are a key permissive factor—Cx45 effectively acting as a brake on contractile activity in the uterus. As parturition approaches, Cx45 becomes downregulated, allowing contractions to increase in strength and frequency.

## Acknowledgements

R.E.S. is grateful to EPSRC for a PhD studentship through the MOAC Doctoral Training Centre at the University of Warwick: grant number EP/F500378/1. C.M. gratefully acknowledges funding for a PhD studentship from the BBSRC through the Systems Biology Doctoral Training Centre at the University of Warwick: grant number BB/G530233/1.

## Appendix A. Derivation of the FitzHugh–Nagumo model

We follow the derivation provided by Keener & Sneyd [27]. Consider the circuit given in the electronic supplementary material, figure S2*b* showing a simplified model of the cell membrane. From Kirchoff's laws, one obtains
A1aand
A1bwhere *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 electronic supplementary material, figure S2

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

*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 [15]) 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) can then be rewritten as follows
A2aand
A2bwhere , *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
A3

- Received July 4, 2014.
- Accepted September 4, 2014.

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