Spatial Noise in Coupling Strength and Natural Frequency within a Pacemaker Network; Consequences for Development of Intestinal Motor Patterns According to a Weakly Coupled Phase Oscillator Model

Pacemaker activities generated by networks of interstitial cells of Cajal (ICC), in conjunction with the enteric nervous system, orchestrate most motor patterns in the gastrointestinal tract. It was our objective to understand the role of network features of ICC associated with the myenteric plexus (ICC-MP) in the shaping of motor patterns of the small intestine. To that end, a model of weakly coupled oscillators (oscillators influence each other's phase but not amplitude) was created with most parameters derived from experimental data. The ICC network is a uniform two dimensional network coupled by gap junctions. All ICC generate pacemaker (slow wave) activity with a frequency gradient in mice from 50/min at the proximal end of the intestine to 40/min at the distal end. Key features of motor patterns, directly related to the underlying pacemaker activity, are frequency steps and dislocations. These were accurately mimicked by reduction of coupling strength at a point in the chain of oscillators. When coupling strength was expressed as a product of gap junction density and conductance, and gap junction density was varied randomly along the chain (i.e., spatial noise) with a long-tailed distribution, plateau steps occurred at pointsof low density. As gap junction conductance was decreased, the number of plateaus increased, mimicking the effect of the gap junction inhibitor carbenoxolone. When spatial noise was added to the natural interval gradient, as gap junction conductance decreased, the number of plateaus increased as before but in addition the phase waves frequently changed direction of apparent propagation, again mimicking the effect of carbenoxolone. In summary, key features of the motor patterns that are governed by pacemaker activity may be a direct consequence of biological noise, specifically spatial noise in gap junction coupling and pacemaker frequency.

Introduction

A few consistent observations have been made of slow waves, or the contractions they drive, in the small intestine. Commonly they travel in an aboral direction; their frequency decreases aborally; this frequency gradient often consists of a series of flat plateaus bounded by abrupt steps; clamping the intestine, so as to block conduction, decreases the frequency aboral to the clamp (e.g., Alvarez, 1914; Hasselbrack and Thomas, 1961; Christensen et al., 1966; Diamant and Bortoff, 1969; Szurszewski et al., 1970). The first of these has an obvious and critical function, to encourage food to travel aborally rather than orally along the gut. Early on it was speculated that the other three are aspects of the mechanism underlying the first (Alvarez, 1928). This mechanism is now understood as coupled oscillator theory (Pavlidis, 1973; Winfree, 1980).

Slow waves are generated in the small intestine by a network of interstitial cells of Cajal (ICC-MP) that covers the entire length and circumference of the intestine at the level of the myenteric plexus (Sanders et al., 2014) with a density of around 1000 cells per mm² (Mei et al., 2009). Each ICC-MP can generate a rhythmic depolarization in isolation (Koh et al., 1998; Thomsen et al., 1998). In coupled oscillator theory each element of a network oscillates in isolation but this oscillation is influenced by the oscillations of its neighboring, coupled oscillators. In general two coupled oscillators will tend to minimize both their phase difference, a process called synchronization, and their frequency difference, called entrainment. Synchronization was first observed by Huygens in the seventeenth century between pendulums (Bennett et al., 2002) and explains how networks of oscillators can generate waves. The wave does not propagate, instead it is a coordinated phase difference across the oscillator network, a “phase wave.” When all the oscillators have the same frequency, the phase difference can go to zero (complete synchronization) and the phase wave will have an infinite velocity. If there is a gradient in frequency across the network then the higher frequency oscillators will always be ahead of the phase of the lower frequency oscillators and so the phase wave will appear to propagate from the high to the low frequency part of the network. The velocity of the phase wave will vary, as the phase difference between different frequency oscillators necessarily changes with time. The frequency difference can be minimized by entrainment. Here a distinction must be made between the oscillator's natural frequency, the frequency at which it oscillates in isolation, and the frequency it oscillates at in a network. Entrainment is when an oscillator with a higher natural frequency, “pulls” up the oscillation frequency of an oscillator with a lower natural frequency, toward its own frequency. How much depends on the strength of coupling and difference in natural frequency between the oscillators. In a chain of oscillators with a natural frequency gradient, entrainment pulls a length of the chain to the same frequency, a frequency plateau. A series of plateaus can result, with intervening steps where the frequency difference is too great for entrainment. At a frequency step one of every few phase waves will terminate as a “dislocation.” For instance if plateau A and its neighbor B have wave intervals of 1.5 and 1.4 s, respectively, then every 15th wave in A will terminate at the step between A and B (Parsons and Huizinga, 2015a).

Coupled oscillator theory coalesced from many different components. Of these, van der Pol's work on relaxation oscillators was particularly important. In a relaxation oscillation the amplitude increases gradually for part of its cycle (relaxes) and then abruptly increases before decreasing again. Van der Pol discovered entrainment in electrical circuits that acted as relaxation oscillators and described a simple differential equation that could model this (van der Pol, 1926; van der Pol and van der Mark, 1927). He also noted that the rhythmicity of the heart had the properties of a relaxation oscillator (van de Pol, 1940). In 1961 Fitzhugh showed that the Hodgkin Huxley equations could be simplified to a second order differential equation that contains the van de Pol equation as a limiting case, now called the Fitzhugh-Nagumo or Bonhoeffer-van der Pol equation (Fitzhugh, 1961). In the same year Bortoff introduced coupled oscillator theory to explain the properties of slow waves in the small intestine (Bortoff, 1961). A decade later Diamant showed that a chain of Fitzhugh-Nagumo oscillators could reproduce the frequency plateaus of the small intestine (Diamant et al., 1970). This was followed by many studies on chains of coupled relaxation oscillators by gastrointestinal physiologists (Brown et al., 1971, 1975; Sarna et al., 1971, 1972; Specht and Bortoff, 1972; Robertson-Dunn and Linkens, 1974; Akwari et al., 1975; Linkens, 1977; Linkens and Datardina, 1977; Publicover and Sanders, 1989; Daniel et al., 1994; Aliev et al., 2000; Gizzi et al., 2010). In 1967 Winfree published his insight that when coupling between oscillators is “weak,” in the sense that oscillators only influence each other's phase and not their amplitude (no oscillator deviates significantly from its limit cycle), then the dynamics of any such oscillator can be described by a single differential equation of the oscillator's phase (Winfree, 1967, 1980). This idea was developed over the next two decades into what is now known as the theory of weakly coupled oscillators (Schwemmer and Lewis, 2012). It provides a practical means to “phase reduce” any model, biophysical or otherwise, of two or more coupled oscillators to a single first order differential equation (Equation 1). This equation has only two variables, natural frequency and coupling strength, and a function, the interaction function. Significantly for the physiologist all but one of these can be determined from experimental data. Natural frequency can be measured from the uncoupled oscillation frequency or estimated from the coupled gradient of oscillation frequency. For the interaction function, a stimulus, such as a current pulse for electrical oscillators, is given at different phases of the oscillator's cycle and the change in phase induced by this stimulus measured (Figure 1A). This gives a phase response curve which can be used to derive the interaction function (Figures 1B,C).

FIGURE 1

Figure 1. Phase response curves and the interaction function. (A) An oscillator is represented by a horizontal time line and “fires” at each solid black uptick. The firing time could be defined as the upstroke of the slow wave. The oscillator has a natural period of T. T and all other variables can be defined either in units of radians (T = 2π), seconds (T = some value) or normalized to T (T = 1). A stimulus is given ϕ after a fire and the next fire is τ after this fire. Further measurements can be made, grouped according to two standpoints: (1) Upper time line, the oscillator “should have” fired T after the stimulus (hollow uptick). ϕ′ is the phase difference between when the oscillator did fire and when it “should have.” If ϕ′ = T, the oscillator fired immediately upon stimulation—the phase of the oscillator was fully advanced (box, center of B). Conversely if ϕ′ = 0, the oscillator fired at the time it “should have”—the phase of the oscillator was fully delayed (box, center of B). If ϕ′ = ϕ, then the stimulus had no effect on the firing of the oscillator—it was refractory. Another measurement used is the cophase (θ) the time between the stimulus and the first fire. (2) Lower time line, the oscillator “should have” fired T after it's last fire (hollow uptick). ϕ* is the phase difference between when the oscillator did fire and when it “should have.” (B) There are a number of conventions for the plotting of phase response curves. The abscissa is always the phase of the stimulus (ϕ), but the ordinate can be ϕ′, −ϕ*, ϕ*, or τ (anticlockwise from left). ϕ′(ϕ) is known as an “old phase—new phase” phase response curve (Winfree, 1980). All the phase response curves can be related by simple geometric transforms (arrows). The interaction function, H(ϕ) is a reversed and smoothed version of the infinitesimal phase response curve −ϕ*(ϕ) (Equation 5). (C) τ(ϕ) with T-normalized units for slow waves recorded from the rabbit small intestine, modified from Figure 3A of Cheung and Daniel (1980). Triangles and circles are circular and longitudinal muscle responses, respectively, to 5 ms, 80 V pulses. The slow waves are refractory until halfway through their cycle at which point they phase advance, reaching full phase advance (responding immediately to a stimulus) about 7/10 of the way through the cycle. The blue dashed line is the infinitesimal phase response curve used for our model (Equation 6). (D) Model infinitesimal phase response curve Z(t), impulse V(t) and interaction function H(ϕ_ij). a, b, and h are the parameters used to define Z(t) (Equation 6).

In a previous paper we studied contraction waves in the murine small intestine using the technique of diameter mapping (Parsons and Huizinga, 2015a). We observed frequency plateaus and dislocations, and showed how these were affected when coupling was reduced by the gap junction blocker carbenoxolone. Our aim here was to model this data using a chain of weakly coupled phase oscillators. All but one parameter of the model were fixed using experimental physiological data. We discovered that spatial noise in coupling strength and natural frequency of the pacemaker network play a decisive role in the orchestration of intestinal motor patterns. Part of this work has been published previously in abstract form (Parsons and Huizinga, 2015b).

Materials and Methods

Experiments and Analysis

All procedures were approved and carried out in accordance with regulations of the Animal Ethics Board of McMaster University. Organ bath video recordings of the murine small intestine were made and diameter maps (DMaps) created, as described previously (Parsons and Huizinga, 2015a). A DMap is an image where the x axis corresponds to time, the y axis to distance along the length of the intestine (in all maps shown, top to bottom corresponds to proximal to distal). The image intensity at {x, y} corresponds to the intestine diameter at that time and distance, with black corresponding to contracted and white, relaxed.

The tone of the small intestine usually varies along its length and so there are variations in DMap intensity along its y/distance axis. To enhance the contrast of contractions and aid comparison with DMaps simulated by the model, intensity was equalized across the y axis. For each point on the y axis the mean intensity was calculated across x and this was subtracted from all pixels at that point.

To measure the frequency of waves in DMaps, experimental or modeled, autocorrelation was calculated for each point on the y/distance axis for lags up to 10 s. This gave a {lag, y} image with intensity on a blue-white-red scale, where blue corresponds to negative correlation, white near to zero and red positive correlation. Autocorrelation is always +1 at zero lag, then proceeds to oscillate within the bounds of −1 and +1. The first positive peak, after zero-lag, corresponds to the correlation between every consecutive wave, the second peak to the correlation between every other wave, etc. The lag of each peak is the interval between waves.

Model

Different authors in the weakly coupled oscillator literature use different terminology, sign conventions and symbols for the same equations. In the following we stay as close as possible to the usage of Schwemmer and Lewis (2012). The model consisted of a chain of n coupled phase oscillators. The rate of change of phase (θ) of the ith oscillator is,

$\begin{array}{l}\frac{d{\theta }_i}{dt}={\omega }_i+\frac{1}{2}{\displaystyle \sum _j}{k}_{ij}H\left({\varphi }_{ij}\right)\\ i=1\dots .n\\ {\omega }_i=2\pi ∕{\zeta }_i\\ j=i-1,i+1\\ {\varphi }_{ij}=\left({\theta }_j-{\theta }_i\right)mod2\pi (1)\end{array}$

where ω is the natural frequency (radians/s); ζ is the natural interval (s); k is the coupling strength (unitless); H is the interaction function (radians/s). n was 200, large enough to give a fairly fine grained simulation and small enough for a fast numerical solution.

The model was solved in MATLAB (MathWorks, Natick, MA) using the ode45 Runge-Kutta method. It was found that the solution was sensitive to rounding errors in θ as θ became larger. To avoid these errors the model was solved in consecutive blocks of 100/max(ω) seconds. The 2π modulo of the last phases of one block were used as the initial phases for the next block (the initial phases for the first block were random). In this way θ never exceeded much beyond 100. With a Intel Core i5-2300 (2.8 GHz) processor it typically took 3–4 s for a 250 s simulation of 200 oscillators. Solutions were returned at an interval of 33 ms, the video sample interval for the experimental DMaps. DMaps (with a y/length axis of n pixels) were created from the solution by converting phases to amplitudes by sine function.

The natural interval (ζ) was normally distributed with standard deviation σ_ζ about a linear gradient along the chain from ζ_low to ζ_high,

$\begin{array}{l}{\zeta }_i={\zeta }_{low}+i\left({\zeta }_{high}-{\zeta }_{low}\right)∕n+{\sigma }_{\zeta }U(2)\end{array}$

where U is a standard normal variate. Based on data from the murine small intestine, ζ_low was 1.2 s and ζ_high was 1.5 s (Figure 4A of Parsons and Huizinga, 2015a). σ_ζ was varied (see Section Results).

Coupling strengths (k) were symmetrical (k_ij = k_ji). In the simplest model, k was the same value for all oscillators except at selected points (Figures 2–4). In a more complex model (Figures 6, 7), k was the product of a gap junction conductance (g_junc) and gap junction density (d_junc). A decrease in g_junc modeled gap junction blockers. d_junc had a reversed Lévy distribution restricted between 0 and 1,

$\begin{array}{l}{k}_{ij}=k_{ji}=g_{junc}{d}_{junc}\\ d_{\text{junc}}=1-X(3)\end{array}$

where X was a random variate from the Lévy stable distribution (Weron and Weron, 2005) with the characteristic function,

$\begin{array}{l}\phi \left(t\right)=exp\left[-{\sigma }^{\alpha }|t{|}^{\alpha }\left(1-i\beta \text{sgn}(t)tan(\pi \alpha ∕2)\right)+it\mu \right](4)\end{array}$

with characteristic exponent α = 0.5, skewness β = 0.5, scale s = 0.02, and location μ = 0. Random variates were generated with the code of Veillette (2008) until a variate was in the range (0, 1).

FIGURE 2

Figure 2. Phase wave patterns and frequency as coupling strength (k) is varied. Model diameter maps (left), and their autocorrelation (right). The top of the maps correspond to the proximal (higher natural frequency) end. Time scales are indicated below (C). Autocorrelation lag runs from left to right and red is positive correlation while blue is negative. After zero lag the first peak in correlation (first red band) corresponds to correlation between every consecutive oscillation, the second peak (second red band) corresponds to correlation between every other oscillation and so on (see Section Materials and Methods). The model has a noiseless frequency gradient (ζ_low = 1.2 s, ζ_high = 1.5 s, σ_ζ = 0 s) and uniform coupling strength. (A) k = 3. (B) k = 1. (C) k = 0.1.

FIGURE 3

Figure 3. Induction of frequency steps by localized decoupling. (A) Model diameter map (left), and its autocorrelation (right). The top of the map corresponds to the proximal (higher frequency) end. Time scales are indicated below (B). The model has the same parameters as in Figure 2A, but with coupling strength reduced by 90% (k = 0.3) at the points marked with red arrows. At the points of reduced coupling there are clear steps in frequency, corresponding to positions of dislocations in the DMap. (B) Experimental diameter map of murine small intestine (left) in the presence of 0.5 mM lidocaine to block enteric nervous system activity and its autocorrelation (right). The top of the map corresponds to the proximal end. (C) A dislocation from the model DMap in (A). (D) A dislocation from the small intestine DMap in (B).

FIGURE 4

Figure 4. Dependence of dislocations on strength of localized decoupling. All 200 oscillators had the same coupling strength (k = 3), except between oscillators 100 and 101 where the coupling was reduced to a percentage of this, indicated at left. Left, DMaps showing the middle 20 oscillators. Red ticks indicate times of wave drops (dislocations). Right, corresponding autocorrelations.

The interaction function (H) is periodic over the period T = 2π radians. It is proportional to the one sided convolution over T, of the infinitesimal phase response curve of the oscillator (Z) and the impulse of the neighboring oscillator (V),

$\begin{array}{l}H\left({\varphi }_{ij}\right)=\frac{1}{T}{\displaystyle \underset{0}{\overset{T}{\int }}Z(t)\left[V(t+{\varphi }_{ij})-V(t)\right]dt}\\ =\frac{1}{T}{\displaystyle \underset{0}{\overset{T}{\int }}Z(t)V(t+{\varphi }_{ij})dt}-\frac{1}{T}{\displaystyle \underset{0}{\overset{T}{\int }}Z(t)V(t)dt}\\ =\frac{1}{T}{\displaystyle \underset{0}{\overset{T}{\int }}Z(t)V(t+{\varphi }_{ij})dt}-C\\ H\left(-{\varphi }_{ij}\right)=\frac{1}{T}{\displaystyle \underset{0}{\overset{T}{\int }}Z(t)V(t-{\varphi }_{ij})dt}-C(5)\end{array}$

Thus, if the impulse is a square wave, the interaction function is a reversed and smoothed version of the infinitesimal phase response curve, with the impulse as the “boxcar” smoothing kernel (Figure 1D). When the impulse is a delta function, H(−ϕ_ij) ~ Z(ϕ_ij).

We based the model's infinitesimal phase response curve on the experimental phase response curve of slow waves in the rabbit small intestine determined by Cheung and Daniel (1980) (Figure 1C). Ideally the infinitesimal phase response curve represents the phase response to an infinitely short stimulus (a delta function). The stimulus used by Cheung and Daniel was 5 ms, 1/740 of the 3.7 s average cycle length, so we felt brief enough for approximation. The infinitesimal phase response curve was,

$\begin{array}{l}Z(t)=\{\begin{array}{cc}0& t<a\\ h\frac{t-a}{b-a}& a<t<b\\ -h\frac{x-2\pi }{2\pi -b}& t>b\end{array}\\ a=\pi \\ b=1.3\pi \\ h=2\pi -b(6)\end{array}$

where a is the length of the cycle when the oscillator is refractory to a stimulus and; b is the point in the cycle beyond which the oscillator is fully phase advanced by a stimulus (i.e., responds to a stimulus by immediate “firing” or fully synchronizes with the stimulating oscillator) (Figure 1D).

Physiologically the impulse is the part of the ICC depolarization that can effect a neighboring, coupled ICC (it doesn't have to correspond to the whole slow wave, just up to where there is no more phase change induced). This was modeled as a simple square wave of width w, as a fraction of the complete cycle,

$\begin{array}{l}V(t)=\{\begin{array}{cc}{(2\pi w)}^{-1}& t\le 2\pi w\\ 0& t>2\pi w\end{array}\\ w=0.1(7)\end{array}$

Z(t) and V(t) were one side convolved to create a “look up table” of H(ϕ_ij) at discrete values of ϕ. This was passed to the ODE evaluation function and H for a particular value of ϕ was found from it by linear interpolation.

Results

Frequency Steps and Localized Decoupling

The parameters of the model were determined as much as possible by experimental data. The natural frequency gradient was estimated from the observed contraction intervals at either end of the mouse intestine, 1.2–1.5 s, proximal to distal (Parsons and Huizinga, 2015a). The interaction function was calculated from the phase response curve of slow waves in the rabbit small intestine (Cheung and Daniel, 1980) and a square impulse representing the depolarising phase of the slow wave. This left one free parameter, the coupling strength, k.

When k was varied over a wide range, phase waves propagated distally along the length of the chain but there were no frequency steps or dislocations (Figure 2). At higher coupling strengths the whole chain was entrained to the same frequency (Figure 2A). The lack of frequency steps was not due to coupling being so strong as to result in entrainment to a single oscillator (the proximal, highest frequency oscillator). At lower coupling strengths entrainment weakened but instead of steps and plateaus forming, frequency varied smoothly, eventually following the natural frequency gradient (Figures 2B,C). Also the phase wave velocity increasingly oscillated as coupling was reduced. Smooth variations in wave frequency or velocity of this magnitude were not seen in the mouse small intestine (Parsons and Huizinga, 2015a and Figure 3B). In previous coupled oscillator models frequency steps occurred spontaneously given a large enough frequency gradient (Diamant et al., 1970; Brown et al., 1971, 1975; Sarna et al., 1971, 1972; Specht and Bortoff, 1972; Robertson-Dunn and Linkens, 1974; Akwari et al., 1975; Linkens, 1977; Linkens and Datardina, 1977; Publicover and Sanders, 1989; Kopell et al., 1990; Daniel et al., 1994). It appeared that the physiological frequency gradient in the mouse was not large enough to explain frequency steps. Instead we made the novel hypothesis that steps in the mouse result from localized decoupling between ICC at the frequency step. Reducing coupling between two oscillators would allow a distal (lower natural frequency) oscillator to escape entrainment by its proximal neighbor and thus allow a plateau step to form. Indeed when this was done by reducing k by 90% at two points, two frequency steps (three plateaus) formed and dislocations occurred rhythmically at these points (Figure 3A). This pattern was very similar to that seen in diameter maps of the whole mouse small intestine in vitro (Figure 3B). Wave velocity decreased distally across each plateau in both the model and intestine (Figures 3A,B) as phase lag increased due to the larger difference between the oscillator's natural and entrained frequency (Somers and Kopell, 1995). Dislocations were also very similar between the intestine and model (Figures 3C,D).

As k was incrementally reduced between the middle two oscillators of the chain, there was no distinguishable change until k was reduced to 40% of the control value of three, when there was an instantaneous change in phase at the point of reduction, a “phase slip.” (Figure 4). As k was reduced further the magnitude (phase difference) of this slip increased, until at 20% the waves proximal and distal to the reduction point were anti phase. At 17.5% a frequency step and wave drops were induced. With further reduction down to 0%, the frequency of wave drops increased as the magnitude of the step increased.

Spatial Noise in Coupling Strength

In our previous study of the small intestine, frequency steps were very stable in position. Inhibition of gap junctions with carbenoxolone induced new steps but the old steps were conserved (Figures 6B,E in Parsons and Huizinga, 2015a). This suggests to us that step positions reflect underlying structural discontinuities rather than being the result of a physiochemical dynamic equilibrium. The model suggests that these discontinuities are a reduction in coupling strength (Figure 3), either in the density of gap junctions or equivalently the density of ICC or cellular connections between ICC. It is possible that discontinuities occur at well spaced points against a uniform background as modeled above (Figure 3A), perhaps through a targeted developmental mechanism. Or perhaps more likely the coupling strength could vary randomly along the length of the intestine, at certain points falling below the threshold to give steps. In other words the coupling strength would be spatially noisy.

If coupling strength is spatially noisy, its distribution must be known to model it. Given the low number of steps under basal conditions and that their number only increases 2–3 times with carbenoxolone (Figure 6 in Parsons and Huizinga, 2015a) we hypothesized that the distribution of coupling strengths should have a long left tail (Figure 5). That way if the threshold for step formation is on the tail, the vast majority of coupling strengths will be above threshold and the relative position of the tail and threshold can vary widely without greatly changing the fraction of strengths below threshold, the number of steps (Figure 5B). We put these considerations into a model for the coupling strength (Equation 3) (Figure 5). k was the product of a gap junction conductance (g_junc) and a gap junction density (d_junc). d_junc was varied randomly along the chain between 0 and 1 according to a Lévy distribution, a long tailed distribution. g_junc was a single value that was decreased to model inhibition of gap junction conductance. Steps occurred at points of low d_junc (Figure 6A). As g_junc was decreased six fold, the number of steps/plateaus increased from ~1–2 to 5–10 (Figures 6B,C), in line with the increase seen with carbenoxolone in the small intestine (Figure 6 in Parsons and Huizinga, 2015a).

FIGURE 5

Figure 5. Coupling strength (k) distribution and the number of frequency steps. The number of steps is proportional to the area of the k distribution left of the threshold (t) for inducing steps. As the distribution of k is contracted P(k) → P(k/a) or shifted left P(k) → P(k − a), the area left of threshold increases and so does the number of steps. Starting with normal (A) and long left tailed (B) distributions with equal area left of threshold (dark shading), when each is contracted by the same amount the area left of threshold (light shading) is much greater in the normal distribution. The same relative effect occurs with shift. Thus, the parameters of a normal distribution have to be much more finely tuned to obtain a particular number of steps. This is important when considering coupling strength in the small intestine because the number of steps does not vary widely between intestines, nor does the number of steps increase exponentially with reduction of coupling by carbenoxolone. In our model k is the product of a gap junction density (d_junc), which varies between oscillators, and a gap junction conductance (g_junc), a single value for all oscillators. Thus, the distribution of k is an expansion/contraction of the distribution of d_junc, according to the value of g_junc. (C) The Lévy-stable distribution of d_junc used in our model (Section Materials and Methods) has a long left tail (black line) in comparison to a normal distribution (gray line). Both distributions were calculated numerically (n = 10⁵). Based on the data in Figure 3 the threshold for step induction is k = 0.525 (17.5% of g_junc = 3). (D) Lévy-stable and Gaussian d_junc distributions were adjusted so that the corresponding k distributions had equal areas left of threshold [P(k < 0.525)] at g_junc = 3, so starting with the same number of steps [N(k < t)] (two, a typical number without carbenoxolone). As g_junc is decreased from three and the distribution of k is contracted, N(k < t) rises much faster for the normal distribution (gray line) than for the Lévy-stable distribution (black line). Both reach P(k < t) = 1, as g_junc goes to t (the whole of the k distribution is to the left of t). Similarly as g_junc is increased above three, N(k < t) decreases further for the normal distribution, but only decreases a little for the Lévy-stable distribution. Again the conclusion is that g_junc has to much more finely tuned with a normal distribution of d_junc then with a long-left tailed distribution.

FIGURE 6

Figure 6. Increase in frequency steps as gap junction conductance is decreased. Model d_junc (left), DMaps (middle), and autocorrelation (right). The top of the maps correspond to the proximal (higher frequency) end. Gap junction density (d_junc) had a long tailed distribution and gap junction conductance (g_junc) was either 3 (A), 1 (B), or 0.5 (C). Red and black horizontal bars, frequency step positions.

Spatial Noise in Natural Frequency

Carbenoxolone produced a zigzagging reversal of wave direction (Figure 7D here and Figure 6A in Parsons and Huizinga, 2015a), but this was not observed when we decreased g_junc (Figure 6). In the presence of carbenoxolone, “high frequency islands” (plateaus with higher frequency than both their distal and proximal neighbors) were seen (Parsons and Huizinga, 2015a) and this gave us a hypothesis as to the development of the zigzagging waves. We hypothesized that natural frequency is spatially noisy, but the magnitude of this noise is normally insignificant compared to the natural frequency change over the span of a plateau. As plateaus get shorter with gap junction inhibition, their natural frequency span also gets smaller, to within the range of the natural frequency noise. Therefore, the chances increase that a local spike in natural frequency will be the highest frequency within its plateau, entraining the slow wave frequency within its plateau, and be higher than any natural frequency in the neighboring proximal plateau, resulting in a high frequency island. Spatial noise in natural frequency would also explain the zigzagging seen with carbenoxolone as local gradients in natural frequency coupled with local variations in coupling strength drive the phase wave in different directions. To test our hypothesis, the natural intervals were varied randomly with a normal distribution of 30 ms standard deviation, about the 1.2–1.5 s proximal to distal gradient. The magnitude of this variation is comparable to the temporal variation of slow wave frequency measured experimentally. It gives a coefficient of variation (c_v = standard deviation/mean) of 0.022 (0.03/1.35). In the dog small intestine (Table 1 of Szurszewski et al., 1970) slow wave frequency for a single animal was 14.8 ± 0.08 S.E. min⁻¹ (measured as slow waves over a minute period), n = 17 and so c_v = 0.08√17∕14.8 = 0.022. As g_junc was decreased the number of plateaus increased as before, but also the phase waves zigzagged and high frequency islands developed (Figures 7A–C). Waves formed a V across each plateau (Figure 7C) exactly like the small intestine in the presence of carbenoxolone (Figure 7D). Also the frequency between steps was often not uniform, varying smoothly in the manner of an oscillating slope rather than a plateau (Figure 7C) and again this was seen in the small intestine (Figure 7D). In some cases frequency steps did not correspond with minima in junction density (Figure 7C).

FIGURE 7

Figure 7. Effect of natural frequency spatial noise on emergence of frequency steps. (A–C) Model d_junc (left), DMaps (middle) and autocorrelation (right). The top of the maps correspond to the proximal (higher frequency) end. The natural intervals were distributed normally about their gradient (1.2–1.5 s) with a standard deviation of 30 ms. Gap junction density (d_junc) had a long tailed distribution and gap junction conductance (g_junc) was either 3 (A), 1 (B) or 0.5 (C). Red and black horizontal bars, frequency step positions. Bars with asterisk indicate where frequency steps are not coincident with dips in d_junc. High frequency islands are indicated by vertical bars at the right of the autocorrelation. (D) DMap and autocorrelation of murine small intestine in the presence of 40 μM carbenoxolone. The top of the map corresponds to the proximal end.

Discussion

Random variation, spatial, or temporal noise, is a natural property of any biological system. The role of temporal noise in slow wave generation has been recognized in the form of stochastic miniature depolarizations called unitary events or spontaneous transient depolarizations (Hashitani et al., 1996; Edwards et al., 1999; Hirst and Edwards, 2001). Here we show that spatial noise in key parameters of the pacemaker network play a decisive role in the orchestration of intestinal motor patterns. The evidence was provided by using experimental data to constrain the few parameters of a chain of weakly coupled phase oscillators and introducing spatial noise.

There are many mechanisms through which a particular distribution of coupling strength could be achieved. Long tailed distributions of connectivity (node degree) in mathematical networks (graphs) can be produced by targeted removal or addition of connections (edges) at nodes that already have few or more connections, respectively (Barabasi and Albert, 1999; Barrat et al., 2008), i.e., positive feedback. Targeted addition is called preferential attachment. A long left tailed distribution of coupling strength could result if connections were removed between already poorly connected ICC (low d_junc). In the brain an excess of neural connections are formed during embryonic development and this is followed postnatally by targeted removal of synapses, dendrites, and neurons themselves, according to each neuron's activity (Hua and Smith, 2004) and thus dependent on the number of connections. This was first observed by Cajal (Schuldiner and Yaron, 2015). He called it “process resorption” and it is now called “pruning.” Gao and colleagues have shown that pruning occurs in the myenteric ICC network (Gao et al., 2013a, 2014). ICC density reaches a peak 2 weeks after birth and after this the distance between ICC increases. Though it is not known whether this pruning is targeted, the fact that it occurs 2 weeks after birth, by which time ICC generate slow waves, suggests the possibility of activity-dependent preferential detachment. It is also possible that gap junction density spatial noise could be independent of a variation in ICC density, ICC having varying gap junction expression levels. There is no data on this. Or gap junction density may not itself be noisy but rather the gap junction conductance. Conductance is determined by the unitary conductance and the open probability and these are controlled by pH and connexin phosphorylation (Nielsen et al., 2012) which could remain stable over minutes.

Imtiaz et al. (2006) applied weakly coupled oscillator theory to slow waves starting with a biophysical model of a single oscillator based on a decade of experimental dissection of the slow wave by van Helden, Hirst, Kito and Suzuki (Hirst and Ward, 2003; van Helden et al., 2010 for review). Each oscillator has two pools of calcium, the cytosol and an intracellular store. Release from the store to the cytosol is through inositol trisphosphate (IP₃) gated channels, with the synthesis of IP₃ modulated by membrane potential. A rise in cytosolic calcium depolarises membrane potential (such as through calcium activated chloride channels) and causes further store calcium release (i.e., calcium induced calcium release). In this way oscillations are generated by feedback between calcium, IP₃ and membrane potential. Imtiaz et al. determined the −ϕ*(ϕ) phase response curve for a single oscillator stimulated with a depolarising pulse (Figures 2, 14 of Imtiaz et al., 2006) and this was almost identical to the experimental phase response curve of Cheung and Daniel (1980), used here. The model cell was refractory for the first half of its cycle, dipped to a small phase delay near T/2, after which the cell rapidly switched to full phase advance. When two oscillators were coupled they synchronized (Imtiaz et al., 2006) and when coupling was lowered in the middle of a chain of oscillators, the oscillators at either side would desynchronize (van Helden and Imtiaz, 2003). This modeled an ingenious experiment where they recorded slow waves from two ends of a strip of gastric muscle and then applied glycyrrhetinic acid (an analog of carbenoxolone) or 2-APB (an IP₃R inhibitor) to the intervening portion as a perfused stream (van Helden and Imtiaz, 2003). In the model the sensitivity of stores to IP₃ was varied randomly across the array of cells with a normal distribution. As this sensitivity largely determined natural frequency it parallels, and could be a biological basis for, the spatial noise of the natural frequency in our model.

Imtiaz et al. also analyzed the interaction function by phase reducing their biophysical model (Imtiaz et al., 2006). The rate of change of the phase difference (ϕ) between two symmetrically coupled, equal frequency oscillators as a function of ϕ is

$\begin{array}{l}\frac{d\varphi }{dt}=\frac{d{\theta }_i}{dt}-\frac{d{\theta }_j}{dt};\varphi ={\theta }_i-{\theta }_j\\ =k\left[H\left({\theta }_j-{\theta }_i\right)-H\left({\theta }_i-{\theta }_j\right)\right]\\ =k\left[H(-\varphi )-H(\varphi )\right]\\ =-k{H}_{odd}(\varphi );H_{odd}=H(\varphi )-H(-\varphi )\\ =kG(\varphi );G(\varphi )=-H_{odd}(\varphi )(8)\end{array}$

The odd part of the interaction function (H_odd), or the equivalent growth function (G), give the equilibrium values of ϕ, where G or H_odd is zero with a slope < 0 (G) or > 0 (H_odd). Thus, one can determine not just that the oscillators can synchronize, but also whether they might stably oscillate out of phase in “phase locked” states. Imtiaz et al. found that as the basal rate of IP₃ synthesis increased, H_odd presented a series of phase locked states (Imtiaz et al., 2006). We did not present G(ϕ) for our interaction function because we were primarily concerned with a qualitative analysis and the extra complications of our model (spatial noise and gradients over multiple oscillators) would confound a simple interpretation of G(ϕ). In our model, neighboring oscillators are always incrementally out of phase (the phase waves have a finite velocity) because of the natural frequency gradient. G(ϕ) is not responsible for this, but if it were there would have to be one equilibrium point very close to ϕ = 0. Stable phase differences of larger magnitude that might reflect G(ϕ) equilibrium points at significant distant from ϕ = 0, were seen at phase slips, but these were at points of decoupling, independent of G(ϕ). The magnitude of the phase slip increased incrementally, again independently of any change in G(ϕ), as the transition between complete entrainment (synchrony) and induction of frequency steps (Figure 4).

Ermentrout and Kopell, “partially motivated by certain phenomena observed in mammalian small intestine” (Ermentrout and Kopell, 1984), modeled chains of coupled phase oscillators with a natural frequency gradient (Ermentrout and Kopell, 1984, 1991; Kopell and Ermentrout, 1986, 1990; Kopell et al., 1990). As with earlier papers on chains of coupled relaxation oscillators (see Section Introduction) frequency plateaus arose without decoupling. In these cases the frequency gradient was often in the region of 50% to obtain 2–3 plateaus vs. 20% in our model. This greater frequency difference was probably enough on its own to break entrainment. Somers and Kopell (1995) gave a proof that when a phase oscillator chain is divided into domains that are out of synchrony with each other (such as with frequency plateaus or decoupling) then waves will propagate in both directions from the center of each domain. This phenomenon, which they called fractured waves, can be clearly seen in both our model and the experimental data (Figure 7). For each plateau, with boundaries determined by decoupling, waves spread in a V from the center of the plateau. However, we only saw this in the model when we introduced the natural frequency spatial noise. Somers and Kopell make no mention of this in their proof, however they do say that in addition to being independent domains, the boundary oscillators at neighboring domains must be out of phase. This might depend on increasing their natural frequency difference with noise.

A waxing and waning of contraction amplitude occurs at frequency steps in the small intestine (Diamant and Bortoff, 1969; Suzuki et al., 1986; Parsons and Huizinga, 2015a). This likely results from superposition of waves at either side of the step, an effect well characterized in acoustics (Discussion of Parsons and Huizinga, 2015a). A weakly coupled oscillator model can have little to say about this interpretation as it does not include amplitude as a variable. Oscillator phase can be converted to amplitude by a fixed function (as done here by sine function) but as each oscillator must pass through the full phase range, the amplitude range during each cycle will also always be the same. Waxing waning of contraction amplitude can occur in the absence of frequency steps through an unrelated mechanism, modulation of slow wave amplitude by a lower frequency pacemaker activity (“phase-amplitude modulation”) (Huizinga et al., 2014).

Arrays of cellular automata have been used to model slow wave propagation in the small intestine (Lammers et al., 2011; Gao et al., 2013b). Lammers et al. (2011) found that when an increasing percentage of cells were made inactive, propagation first slowed and then began to fail. The percentage of active cells would be analogous to coupling strength in our model and though we did not measure velocity changes in the model, propagation velocity does decrease with decoupling by carbenoxolone in the small intestine (Parsons and Huizinga, 2015a).

Weakly coupled oscillator theory is extremely powerful as it can distill into a single equation, biophysical models, via phase reduction, and experimental data and classical physiological concepts such as refractory period, via the phase response curve. Thus, one can tie together biophysics, mathematical analysis, experiment, and modeling to reveal connections that might not otherwise be apparent and better allowing one to feed the other. We have demonstrated a novel connection between spatial variation or noise and the physiological and experimental phenomena of frequency step stability, V-waves and high frequency islands, a connection that would not have been in any way obvious or trivial without the guide of the theory. This connection suggests further experimental, biophysical and mathematical investigations.

Author Contributions

SP designed and executed the experiments, and wrote the manuscript. JH supervised and contributed to manuscript writing.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This study was supported by Canadian Institutes of Health Research Grant MOP12874 and Natural Sciences and Engineering Research Council Grant 386877 to JH.

References

Akwari, O. E., Kelley, K. A., Steinbach, J. H., and Code, C. F. (1975). Electric pacing of intact and transected canine small intestine and its computer model. Am. J. Physiol. 229, 1188–1197.

PubMed Abstract | Google Scholar

Aliev, R. R., Richards, W., and Wikswo, J. P. (2000). A simple nonlinear model of electrical activity in the intestine. J. Theor. Biol. 204, 21–28. doi: 10.1006/jtbi.2000.1069

PubMed Abstract | CrossRef Full Text | Google Scholar

Alvarez, W. C. (1914). Functional variations in contractions of different parts of the small intestine. Am. J. Physiol. 35, 177–193.

Google Scholar

Alvarez, W. C. (1928). The Mechanics of the Digestive Tract, 2nd Edn. New York, NY: Hoeber.

Barabasi, A. L., and Albert, R. (1999). Emergence of scaling in random networks. Science 286, 509–512. doi: 10.1126/science.286.5439.509

PubMed Abstract | CrossRef Full Text | Google Scholar

Barrat, A., Barthelemy, M., and Vespignani, A. (2008). Dynamical Processes on Complex Networks. Cambridge: Cambridge University Press.

Bennett, M., Schatz, M. F., Rockwood, H., and Wiesenfeld, K. (2002). Huygens's clocks. Proc. R. Soc. A 458, 563–579. doi: 10.1098/rspa.2001.0888

CrossRef Full Text | Google Scholar

Bortoff, A. (1961). Slow potential variations of small intestine. Am. J. Physiol. 201, 203–208.

Google Scholar

Brown, B. H., Duthie, H. L., Horn, A. R., and Smallwood, R. H. (1975). A linked oscillator model of electrical activity of human small intestine. Am. J. Physiol. 229, 384–388.

PubMed Abstract | Google Scholar

Brown, B. H., Ng, K. K., Kwong, K., Duthie, H. L., Whittaker, G. E., and Franks, C. I. (1971). Computer analysis and simulation of human gastroduodenal electrical activity. Med. Biol. Eng. 9, 305–314. doi: 10.1007/BF02474084

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheung, D. W., and Daniel, E. E. (1980). Comparative study of the smooth muscle layers of the rabbit duodenum. J. Physiol. 309, 13–27. doi: 10.1113/jphysiol.1980.sp013490

PubMed Abstract | CrossRef Full Text | Google Scholar

Christensen, J., Schedl, H. P., and Clifton, J. A. (1966). The small intestinal basic electrical rhythm (slow wave) frequency gradient in normal men and in patients with variety of diseases. Gastroenterology 50, 309–315.

PubMed Abstract

Daniel, E. E., Bardakjian, B. L., Huizinga, J. D., and Diamant, N. E. (1994). Relaxation oscillator and core conductor models are needed for understanding of GI electrical activities. Am. J. Physiol. Gastrointest. Liver Physiol. 266, G339–G349.

PubMed Abstract | Google Scholar

Diamant, N. E., and Bortoff, A. (1969). Nature of the intestinal slow-wave frequency gradient. Am. J. Physiol. 216, 301–307.

PubMed Abstract | Google Scholar

Diamant, N. E., Rose, P. K., and Davison, E. J. (1970). Computer simulation of intestinal slow-wave frequency gradient. Am. J. Physiol. 219, 1684–1690.

PubMed Abstract | Google Scholar

Edwards, F. R., Hirst, G. D., and Suzuki, H. (1999). Unitary nature of regenerative potentials recorded from circular smooth muscle of guinea-pig antrum. J. Physiol. 519, 235–250. doi: 10.1111/j.1469-7793.1999.0235o.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Kopell, N. (1984). Frequency plateaus in a chain of weakly coupled oscillators. SIAM J. Math. Anal. 15, 215–237. doi: 10.1137/0515019

CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Kopell, N. (1991). Multiple pulse interactions and averaging in systems of coupled neural oscillators. J. Math. Biol. 29, 195–217. doi: 10.1007/BF00160535

CrossRef Full Text | Google Scholar

Fitzhugh, R. (1961). Impulses and physiological states in theoretical models of nerve membrane. Biophys. J. 1, 445–466. doi: 10.1016/S0006-3495(61)86902-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, J., Du, P., O'Grady, G., Archer, R., Farrugia, G., Gibbons, S. J., et al. (2013a). Numerical metrics for automated quantification of interstitial cell of Cajal network structural properties. J. R. Soc. Interface 10:20130421. doi: 10.1098/rsif.2013.0421

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, J., Du, P., O'Grady, G., Archer, R., Gibbons, S. J., Farrugia, G., et al. (2013b). Cellular automaton model for simulating tissue-specific intestinal electrophysiological activity. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2013, 5537–5540. doi: 10.1109/EMBC.2013.6610804

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, J., Sathar, S., O'Grady, G., Han, J., and Cheng, L. K. (2014). Developmental changes in postnatal murine intestinal interstitial cell of Cajal network structure and function. Ann. Biomed. Eng. 42, 1729–1739. doi: 10.1007/s10439-014-1021-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Gizzi, A., Cherubini, C., Migliori, S., Alloni, R., Portuesi, R., and Filippi, S. (2010). On the electrical intestine turbulence induced by temperature changes. Phys. Biol. 7:16011. doi: 10.1088/1478-3975/7/1/016011

PubMed Abstract | CrossRef Full Text | Google Scholar

Hashitani, H., Van Helden, D. F., and Suzuki, H. (1996). Properties of spontaneous depolarizations in circular smooth muscle cells of rabbit urethra. Br. J. Pharmacol. 118, 1627–1632. doi: 10.1111/j.1476-5381.1996.tb15584.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hasselbrack, R., and Thomas, J. E. (1961). Control of intestinal rhythmic contractions by a duodenal pacemaker. Am. J. Physiol. 201, 955–960.

PubMed Abstract | Google Scholar

Hirst, G. D., and Edwards, F. R. (2001). Generation of slow waves in the antral region of guinea-pig stomach–a stochastic process. J. Physiol. 535, 165–180. doi: 10.1111/j.1469-7793.2001.00165.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirst, G. D., and Ward, S. M. (2003). Interstitial cells: involvement in rhythmicity and neural control of gut smooth muscle. J. Physiol. 550, 337–346. doi: 10.1113/jphysiol.2003.043299

PubMed Abstract | CrossRef Full Text | Google Scholar

Hua, J. Y., and Smith, S. J. (2004). Neural activity and the dynamics of central nervous system development. Nat. Neurosci. 7, 327–332. doi: 10.1038/nn1218

PubMed Abstract | CrossRef Full Text | Google Scholar

Huizinga, J. D., Chen, J. H., Zhu, Y. F., Pawelka, A., McGinn, R. J., Bardakjian, B. L., et al. (2014). The origin of segmentation motor activity in the intestine. Nat. Commun. 5, 3326. doi: 10.1038/ncomms4326

PubMed Abstract | CrossRef Full Text | Google Scholar

Imtiaz, M. S., Katnik, C. P., Smith, D. W., and van Helden, D. F. (2006). Role of voltage-dependent modulation of store Ca2+ release in synchronization of Ca2+ oscillations. Biophys. J. 90, 1–23. doi: 10.1529/biophysj.104.058743

PubMed Abstract | CrossRef Full Text | Google Scholar

Koh, S. D., Sanders, K. M., and Ward, S. M. (1998). Spontaneous electrical rhythmicity in cultured interstitial cells of cajal from the murine small intestine. J. Physiol. 513, 203–213. doi: 10.1111/j.1469-7793.1998.203by.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kopell, N., and Ermentrout, G. B. (1986). Symmetry and phaselocking in chains of weakly coupled oscillators. Commun. Pure Appl. Math. 39, 623–660. doi: 10.1002/cpa.3160390504

CrossRef Full Text | Google Scholar

Kopell, N., and Ermentrout, G. B. (1990). Phase-transitions and other phenomena in chains of coupled oscillators. SIAM J. Appl. Math. 50, 1014–1052. doi: 10.1137/0150062

CrossRef Full Text | Google Scholar

Kopell, N., Zhang, W., and Ermentrout, G. B. (1990). Multiple coupling in chains of oscillators. SIAM J. Math. Anal. 21, 935–953. doi: 10.1137/0521052

CrossRef Full Text | Google Scholar

Lammers, W. J. E. P., Al-Bloushi, H. M., Al-Eisaei, S. A., Al-Dhaheri, F. A., Stephen, B., John, R., et al. (2011). Slow wave propagation and plasticity of interstitial cells of Cajal in the small intestine of diabetic rats. Exp. Physiol. 96, 1039–1048. doi: 10.1113/expphysiol.2011.058941

PubMed Abstract | CrossRef Full Text | Google Scholar

Linkens, D. A. (1977). The stability of entrainment conditions for RLC coupled Van der Pol oscillators used as a model for intestinal electrical rhythms. Bull. Math. Biol. 39, 359–372.

PubMed Abstract | Google Scholar

Linkens, D. A., and Datardina, S. (1977). Frequency entrainment of coupled Hodgkin-Huxley-type oscillators for modeling gastro-intestinal electrical activity. IEEE Trans. Biomed. Eng. 24, 362–365. doi: 10.1109/TBME.1977.326143

PubMed Abstract | CrossRef Full Text | Google Scholar

Mei, F., Zhu, J., Guo, S., Zhou, D. S., Han, J., Yu, B., et al. (2009). An age-dependent proliferation is involved in the postnatal development of interstitial cells of Cajal in the small intestine of mice. Histochem. Cell Biol. 131, 43–53. doi: 10.1007/s00418-008-0515-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Nielsen, M. S., Axelsen, L. N., Sorgen, P. L., Verma, V., Delmar, M., and Holstein-Rathlou, N. H. (2012). Gap Junctions. Compr. Physiol. 2, 1981–2035. doi: 10.1002/cphy.c110051

PubMed Abstract | CrossRef Full Text | Google Scholar

Parsons, S. P., and Huizinga, J. D. (2015a). Effects of gap junction inhibition on contraction waves in the murine small intestine in relation to coupled oscillator theory. Am. J. Physiol. Gastrointest. Liver Physiol. 308, G287–G297. doi: 10.1152/ajpgi.00338.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Parsons, S. P., and Huizinga, J. D. (2015b). How do ICC Network Disturbances Affect Small Intestine Motor Patterns? Canadian Digestive Diseases Week. Banff, AB: Can J Gastroent Hepatol.

Pavlidis, T. (1973). Biological Oscillators: their Mathematical Analysis. New York, NY: Academic Press.

Google Scholar

Publicover, N. G., and Sanders, K. M. (1989). Are relaxation oscillators an appropriate model of gastrointestinal electrical activity? Am. J. Physiol. Gastrointest. Liver Physiol. 256, G265–G274.

PubMed Abstract | Google Scholar

Robertson-Dunn, B., and Linkens, D. A. (1974). A mathematical model of the slow-wave electrical activity of the human small intestine. Med. Biol. Eng. 12, 750–758. doi: 10.1007/BF02477440

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanders, K. M., Ward, S. M., and Koh, S. D. (2014). Interstitial cells: regulators of smooth muscle function. Physiol. Rev. 94, 859–907. doi: 10.1152/physrev.00037.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarna, S. K., Daniel, E. E., and Kingma, Y. J. (1971). Simulation of slow-wave electrical activity of small intestine. Am. J. Physiol. 221, 166–175.

PubMed Abstract | Google Scholar

Sarna, S. K., Daniel, E. E., and Kingma, Y. J. (1972). Effects of partial cuts on gastric electrical control activity and its computer model. Am. J. Physiol. 223, 332–340.

PubMed Abstract | Google Scholar

Schuldiner, O., and Yaron, A. (2015). Mechanisms of developmental neurite pruning. Cell. Mol. Life Sci. 72, 101–119. doi: 10.1007/s00018-014-1729-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwemmer, M. A., and Lewis, T. J. (2012). “The theory of weakly coupled oscillators,” in Phase Response Curves in Neuroscience, eds N. W. Schultheiss, A. A. Prinz, and R. J. Butera (New York, NY: Springer), 518.

Somers, D., and Kopell, N. (1995). Waves and synchrony in networks of oscillators of relaxation and non-relaxation type. Physica D 89, 169–183. doi: 10.1016/0167-2789(95)00198-0

CrossRef Full Text | Google Scholar

Specht, P. C., and Bortoff, A. (1972). Propagation and electrical entrainment of intestinal slow waves. Am. J. Dig. Dis. 17, 311–316. doi: 10.1007/BF02231730

PubMed Abstract | CrossRef Full Text | Google Scholar

Suzuki, N., Prosser, C. L., and DeVos, W. (1986). Waxing and waning of slow waves in intestinal musculature. Am. J. Physiol. 250, G28–G34.

PubMed Abstract | Google Scholar

Szurszewski, J. H., Elveback, L. R., and Code, C. F. (1970). Configuration and frequency gradient of electric slow wave over canine small bowel. Am. J. Physiol. 218, 1468–1473.

PubMed Abstract | Google Scholar

Thomsen, L., Robinson, T. L., Lee, J. C., Farraway, L. A., Hughes, M. J., Andrews, D. W., et al. (1998). Interstitial cells of Cajal generate a rhythmic pacemaker current. Nat. Med. 4, 848–851. doi: 10.1038/nm0798-848

PubMed Abstract | CrossRef Full Text | Google Scholar

van de Pol, B. (1940). Biological rhythms considered as relaxation oscillations. Acta Med. Scand. 108, 76–88.

Google Scholar

van der Pol, B. (1926). On “relaxation oscillations.” Philos. Mag. 2, 978–992. doi: 10.1080/14786442608564127

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Pol, B., and van der Mark, J. (1927). Frequency demultiplication. Nature 120, 363–364. doi: 10.1038/120363a0

CrossRef Full Text | Google Scholar

van Helden, D. F., and Imtiaz, M. S. (2003). Ca2+ phase waves: a basis for cellular pacemaking and long-range synchronicity in the guinea-pig gastric pylorus. J. Physiol. 548, 271–296. doi: 10.1113/jphysiol.2002.033720

PubMed Abstract | CrossRef Full Text | Google Scholar

van Helden, D. F., Laver, D. R., Holdsworth, J., and Imtiaz, M. S. (2010). Generation and propagation of gastric slow waves. Clin. Exp. Pharmacol. Physiol. 37, 516–524. doi: 10.1111/j.1440-1681.2009.05331.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Veillette, M. (2008). Avaiable online at: http://math.bu.edu/people/mveillet/research.html (Accessed June 1, 2015).

Weron, A., and Weron, R. (2005). “Computer simulation of Lévy α-stable variables and processes,” in Chaos—The Interplay Between Stochastic and Deterministic Behaviour. Lecture Notes in Physics, eds P. Garbaczewski, M. Wolf, and A. Weron (Berlin: Springer), 379–392.

PubMed Abstract

Winfree, A. T. (1967). Biological rhythms and the behavior of populations of coupled oscillators. J. Theor. Biol. 16, 15–42. doi: 10.1016/0022-5193(67)90051-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Winfree, A. T. (1980). The Geometry of Biological Time. Heidelberg: Springer-Verlag.

Google Scholar