1 Introduction

Rhythmic movements in animals are controlled by neuronal networks that eventually make motoneurons fire at appropriate times. Locomotor activities like swimming, flying or walking form a special class of rhythmic movements. Simple vertebrates and invertebrates have proved to be useful in understanding the complex central pattern generators (CPG) that are responsible for these locomotor behaviours (Stein et al. 1997). One of the best characterised CPGs is in young tadpoles of the frog Xenopus laevis and underlies swimming. This simple vertebrate has mainly been studied shortly after hatching (at stage 37/38 of Nieuwkoop and Faber 1956). Recent progress in understanding the organization of the tadpole networks controlling swimming make it timely to construct a full length population model of these networks. Firstly, new anatomical information is available on the distribution and axonal projections of interneurons and motoneurons (Yoshida et al. 1998; Roberts et al. 1999; Li et al. 2001, 2006). Secondly, physiological observations have established the functional longitudinal connection patterns of reciprocal inhibitory interneurons (Soffe et al. 2001), revealed more details of the central projections and connections of the specific class of premotor excitatory interneurons that drive swimming (Li et al. 2006), and shown that motoneurons may make central cholinergic synapses with more caudal neurons and local electrical connections with each other (Perrins and Roberts 1995a, b, c).

Most previous models of longitudinal coordination of locomotor waves have used a series of coupled, segmental oscillators (Tunstall et al. 2002; Hill et al. 2003). However, there is little evidence for segmentation within the spinal cord. To study intersegmental coordination in the lamprey, Wadden et al. (1997) therefore used a computer model that was based on evenly distributed neurons along the body axis with no segmental boundaries or coupled oscillators. In the hatchling Xenopus tadpole there is also little evidence for internal spinal cord segmentation but neurons are not evenly distributed (e.g. Li et al. 2004b; Yoshida et al. 1998). Therefore we wanted to use a continuous model, in which the spinal neurons of the CPG are distributed along the R-C axis in a realistic manner with axon lengths corresponding to projection distances found in anatomy and physiology.

Three forms of co-ordination are necessary during swimming in tadpoles (Kahn and Roberts 1982; Kahn et al. 1982): (1) left and right side muscles have to contract alternately to bend the trunk (frequency: 12–25 Hz), (2) waves of contractions have to propagate in a head-to-tail, rostro-caudal (R-C) direction along the body so the animal is propelled forward in the water (R-C delay: 2–5 ms/mm and does not scale with cycle period, Tunstall and Roberts 1991) and (3) motoneuron firing has to be synchronized locally (ventral root burst durations: 5–10 ms) for more efficient and powerful muscle contractions within a myotome. Since the basic principles of left and right side alternations have been established both experimentally (Dale 1985) and in models (Roberts and Tunstall 1990; Sautois et al. 2007), this study focuses on the longitudinal (or intersegmental) and local (or intrasegmental) coordination necessary for swimming. It uses a continuous population model with a length dimension to investigate the tadpole’s spinal CPG for swimming.

The main questions we address are: (a) What features of neuron population and axon and synapse distribution are required for stable, self-sustaining swimming activity (lasting from a few seconds to minutes, Kahn and Roberts 1982) with a head-to-tail progression of activity suitable to drive swimming? (b) What is the role of the proposed central cholinergic motoneuron-to-interneuron feedback synapses? (c) What is the role of electrical coupling between motoneurons and is this homogeneous along the spinal cord? (d) How sensitive are our findings to changes in cellular properties of neurons?

2 Methods

2.1 Model neurons and the four cell test oscillator

For modelling we used GENESIS software version 2.1 run under Red Hat 7.0 and Debian 3.0 Linux on Pentium IV PCs.

The small and electrotonically compact neurons of the tadpole (Wolf et al. 1998) were modelled as single isopotential compartments. Axons are represented by conduction delays of 3.64 ms/mm and at synapses a 0.5 ms synaptic delay was added (Dale and Roberts 1985). Motoneurons and interneurons have the same cellular properties. We modelled the CPG in the tadpole using two different neurons (see the appendix for the lists of parameters used) to check if the main conclusions of our recent study are sensitive to the cellular properties of neurons. The first type of neuron mimicked the membrane properties and spike features recorded from spinal rhythmic neurons by sharp microelectrodes (Soffe 1990) and produced non-overshooting single spikes to current injections over a wide range of current intensities. The properties of the second type neuron were set to have less negative resting membrane potential and to produce overshooting spikes, the features seen by whole-cell patch electrodes in more recent experiments (Li et al. 2004a). The membrane kinetics and the density of voltage-dependent K+ and Na+ channels of the type 1 neuron was adapted to GENESIS software from the model of Roberts and Tunstall 1990. The type 2 neuron was obtained from the first one by modifying the resting membrane potential, and by changing the density of the voltage-dependent sodium and potassium channels to get overshooting spikes (up to +30 mV) in response to current injections.

The ability of these model neurons to produce patterns of swimming-like neuronal activity if connected in appropriate networks was first checked in a simple four neuron oscillator (see Section 3, Fig. 1; Roberts and Tunstall 1990).

Fig. 1
figure 1

The four cell oscillator network used to test model neurons: with reciprocal inhibition between half-centres and feedback excitation within each half-centre. (a) Activity was initiated by a single sensory pathway excitation (SEPSP) given to the left (LSEPSP) and right (RSEPSP) sides of the spinal cord with 20 ms delay. Closed triangles are excitatory and circles are inhibitory synapses. Synapses onto boxes mean connection to each neuron within box. e excitatory neuron, i inhibitory neuron. (b) Activity of the excitatory neurons on the two sides following SEPSP (for properties of neurons and synapses see Appendix)

2.2 Building the population model

In the tadpole’s spinal cord, three types of interneurons and the motoneurons are known to be active during swimming (Roberts 2000; Li et al. 2004b). The excitatory interneurons (eIN) are the ‘descending interneurons’. The eINs always have a descending axon but some also have an ascending axon and these are common in the hindbrain (Li et al. 2006) They supply dual component excitation on the same side which is mainly glutamatergic with a smaller, recently discovered cholinergic component (Li et al. 2004c). Since the time course of this cholinergic component is similar to that of the fast component of the glutamate excitation (Li et al. 2004c) we did not model it as a separate process and we will refer to these synapses as glutamate synapses. This is a simplifying assumption. The inhibitory interneurons (iIN) are the ‘commissural interneurons’ with axons that cross the cord where they branch to ascend and descend on the opposite side (Yoshida et al. 1998). The iINs use glycine as a neurotransmitter. Since ascending inhibition has little effect on swimming (Roberts and Alford 1986; Green and Soffe 1998), and basic features of swimming could be reproduced in a shorter model without ascending inhibitory neurons (Li et al. 2006), we made another simplifying assumption and decided not to include the recurrent inhibitory ascending interneurons in the first population model. Motoneurons (MNs) have descending axons establishing cholinergic synapses on the ipsilateral side (Roberts et al. 1999).

We extended our 4 cell network to a population model of the hindbrain and spinal cord with R-C length dimension of 3.5 mm (0 mm is rostral end of hindbrain). Our aim was to use data from anatomical studies to estimate the distributions of neuronal somata and axon lengths of excitatory and inhibitory interneurons (eIN and iIN) and motoneurons (MN; Fig. 2; Yoshida et. al. 1998; Roberts et. al. 1999; Li et. al. 2001, 2006). We used the simplest possible linear approximations of the longitudinal distributions of neuronal cell bodies [Fig. 2(a)] (see Appendix for equations). Each side of the spinal cord was divided into 100 μm long bins and the number of neurons of each type in each bin was determined by these linear distributions. This means that the number of neurons in a certain bin was deterministic. Positions of neurons within bins were determined using a random number generator, so were probabilistic. Axon lengths were short (<1 mm) relative to the length of the model (3.5 mm) and were either constant (for eINs: descending = 700 μm, ascending = 500 μm) or determined by linear functions of neuron positions (MNs and iINs) [Fig. 2(b)] (see Appendix for equations).

Fig. 2
figure 2

Distributions of neurons and axon lengths based on anatomy. (a) Circles show anatomical estimates of cell body distributions for motoneurons (MN), excitatory descending interneurons (eIN) and reciprocal inhibitory commissural interneurons (iIN; Roberts et al. 1999; Li et al. 2001; Yoshida et al. 1998). Lines show plots of simple functions used to model these distributions and determine neuron numbers at different longitudinal positions in the population model (0 mm = rostral end of hindbrain. Hindbrain ends at 0.8 mm and spinal cord at 3.5 mm). Dashed lines show where this ‘full length model’ was cut to create a ‘reduced length’ spinal model (see later in text). (b) Points show anatomical measurements of descending (solid) and ascending (open) axon lengths for different soma positions of each neuron type (Roberts et al. 1999 and unpublished data of Dr. Roberts’s laboratory). Lines show plots of functions used to determine axon lengths for each neuron type in the population model

2.3 Rules of connections and synapses

Three types of neurons are represented in the population model: motoneurons (MN), excitatory interneurons (eIN) and inhibitory commissural interneurons (iIN). MNs synapse onto more caudal neurons of each type on the same side. EINs synapse onto each neuron type positioned either caudally or rostrally on the same side. Commissural iINs make synapses on the opposite side with all types of neurons positioned both rostrally and caudally. MNs on the same side make local electrical synapses, in an all-to-all manner, within each consecutive 0.15 mm long segment representing the muscle blocks or myotomes. These rules of interconnectivity among spinal neurons are based on experimental evidence (for reviews see Roberts et. al. 1998, 2008, but see Section 4).

Each chemical synaptic connection was made probablistically. If the axon of the a presynaptic neuron was long enough to pass or reach a certain neuron then a connection was made with a certain defined probability (see Appendix). The 3.5 mm population model with the neuron distributions, axon lengths and connections defined above will be referred as ‘full-length model’ in figure legends.

2.4 Statistical analysis

For statistical analysis Microsoft Excel (Microsoft Corp.) and Past (Hammer et al. 2001) softwares were used. For comparisons of means either the two-tailed t test or the Mann–Whitney test was used depending on whether the criteria for using a t test were met. Normality of distributions and equivalence of variances were tested by the Shapiro–Wilk and F tests. Significance level was chosen to be 0.05 in all statistical tests. In comparisons of R-C delays five consecutive cycles were considered at least 0.5 s after the initiation of neural activity. R-C delays were measured by the slopes of linear regression lines fitted to scatterograms of MN firing times vs position on one side of the spinal region of the model for each cycle separately. In quantitative data the means are followed by SDs.

3 Results

3.1 Neurons and basic network oscillator model

To construct a population model, we developed a single compartment neuron with properties matched to the definitions of Soffe (1990) based on sharp electrode recordings. This seemed a reasonable simplification as Wolf et al. (1998) showed virtually no difference in soma potentials generated by synaptic inputs to the soma or to the short dendrites found in tadpole neurons. Using the parameters in the Appendix we constructed model neurons and found that the properties and responses of such a single compartment neuron model were very close to those of our previous three compartment model (Roberts and Tunstall 1990). The model neurons were then connected to form a 4 neuron oscillator network [Fig. 1(a)] with reciprocal inhibition between each side and with positive feedback within each side (half-centre). When given a brief excitation, this network generated an alternating pattern of activity [Fig. 1(b)] very similar to previous models of tadpole swimming (Roberts and Tunstall 1990; Tunstall et. al. 2002). Activity depends on feedback excitation to neurons on the same side and on reciprocal inhibition between the two sides. The feedback excitation sums with the falling phase of the long duration NMDA dependent EPSP, which outlasts the spike. The timed inhibition from the other side evokes a delayed spike by post-inhibitory rebound (Li et al. 2006). Once this simple four neuron model had been shown to have stable rhythmic activity, we constructed a population model of the spinal pattern generator for swimming in the tadpole.

3.2 Establishing the population model

We gave a length dimension to the model, by distributing inhibitory and excitatory interneurons (iIN and eIN) and the motoneurons (MN) along a 3.5 mm long R-C axis. The build up of synaptic connections (Fig. 3) was automatic and based on the anatomically defined distributions of neurons and axon lengths (Appendix).

Fig. 3
figure 3

Diagram of the full-length population model with the network’s synaptic connections. All connections are symmetrical. Axons (solid lines) can ascend, descend, cross the cord and leave the cord to innervate muscles. Closed and open triangles are glutamate and acetylcholine mediated excitatory synapses, formed by excitatory interneurons (eINs) and motoneurons (MNs), closed circles are glycinergic synapses made by inhibitory interneurons (iINs) and resistance symbols are electrical synapses between intrasegmental MNs. If a synapse is made onto a box then the presynaptic neuron makes synapses onto neurons within the box with a specified probability if the axon passes the neuron soma position. MNs within 150 μm long segments are interconnected by electrical synapses. Activity is initiated by sensory excitation (SEPSP) at fixed times. Big open arrows and scale bar show R-C directions and extent of the population model. The “reduced model” from 1 to 2.5 mm is indicated on far right. Note that some types of connections were not used in certain simulations to reveal their possible significance or to simplify the model. See text for details

3.2.1 Establishing synaptic connection strengths

We followed a series of steps to determine the strengths of synaptic connections in the population model. The first of these was to use the simple four neuron model (Fig. 1) to find the levels of total, excitatory and inhibitory conductances needed for each neuron of the network to generate stable activity following sensory input. The next step was to distribute these conductances within the much larger population network. Because each neuron in the population may receive synapses from many more neurons than in the simple model, the unitary synaptic conductances needed to be much smaller. To achieve this, we first allowed the Genesis simulator to build a population network where the probability of synapse formation was set at 1. The maximum possible number of synapses made onto individual neurons could then be counted. We next set the unitary conductances to lie within the ranges estimated by Dale (1995) for synapses onto motoneurons: glutamate EPSP range 0.3–0.9 nS and glycinergic IPSP range 4–12 nS. We then reduced the probabilities that determined the frequency of synaptic contacts between neurons and passing axons until the total synaptic conductance per neuron was similar to that used in the four neuron model for each type of synapse (final probabilities were: 0.3 for synapses from eINs, 0.2 for synapses from iINs, 0.5 for synapses from MNs). To simplify the population model for initial testing, we did not include cholinergic and electrical connections made by MNs (Perrins and Roberts 1995a, b), and eINs had only descending axons (CF Li et al. 2006).

Once the unitary synaptic conductances were assigned (these were constant for each synapse and each type of neuron), the population model was tested in response to brief sensory excitation. Although it was based on essentially the same connectivity, and although synaptic conductances as well as cellular properties for individual neurons were equivalent to the four-cell model, the population model failed to produce sustained neural activity. Neurons gradually stopped firing and finally all neurons became silent in less than 3 s. In terms of coordination, while the model produced alternating activity between the two sides, neurons failed to establish the rostral to caudal firing sequence on the same side which is found during swimming (Tunstall and Roberts 1991). Reliable alternating activity within the normal swimming frequencies (12–25 Hz), with rostro-caudal delay and sustained for longer than 10 s could not be achieved even by modifying the neuron or axon distributions or the probabilities of synapses and their unitary synaptic conductances within the physiological ranges, or by turning on electrical synapses between MNs.

3.2.2 The anatomical basis for gradients in synaptic drive

Our initial population model lacked two key features of a successful swimming model: it failed to give persistent activity, and it did not show a head-to-tail progression of activity. To address the problem of longitudinal coordination, we asked what features of synaptic drive are responsible for the normal R-C sequence of activity during swimming. Physiological experiments have shown that there is a R-C gradient in the strength of synaptic drive to spinal neurons during swimming (Tunstall and Roberts 1994; Zhao et al. 1998). Furthermore, the R-C sequence of neuron firings is present even when acetylcholine synaptic transmission is blocked (Zhao et al. 1998). This suggested that a R-C gradient in excitatory drive from eINs is sufficient to produce the normal head to tail spread of firing. We therefore investigated the longitudinal distribution of excitation by examining factors affecting the longitudinal distribution of eIN axons, on the assumption that axon numbers will be proportional to synapse numbers. In a simple anatomical model, if all eINs had descending axons long enough to reach the full length of the spinal cord, there would be a tail-to-head gradient in axon numbers [Fig. 4(b, c)], and an inappropriate tail-to-head gradient of excitation. If all eINs instead have only 700 μm descending axons (within the normal range found anatomically: Li et al. 2001, 2006), then axon numbers rise to a peak at 1 mm and then decline caudally [Fig. 4(b, c)].

Fig. 4
figure 4

Effects of eIN axon length and projection direction on eIN axon distribution in a full-length, simple anatomical model. (a) distribution of eIN somata [Fig. 2(a)]. (b) sample of 12 eINs to illustrate long (gray) and 700 μm descending axons (black). (c) if eINs have long descending axons the number of axons increases caudally (gray) but with 700 μm descending axons (black) eIN axons peak at 1 mm and decrease caudally. (d) sample of eINs with 500 μm ascending + 700 μm descending axons. (e) Sample of eINs where all have a 700 μm descending axon but 500 μm ascending axons only present for neurons from 0 to 1 mm. (f) With descending and ascending axons eIN axon numbers peak at 1 mm (solid line) but when only neurons from 0 to 1 mm have ascending axons there is a plateau of axon numbers in the hindbrain (dashed line). Hindbrain is shaded pale gray

In the animal, some 50% of excitatory ‘descending interneurons’ in the hindbrain have ascending axons with a mean length of 500 μm (Li et al. 2006). In the spinal cord only 15% have such axons and their lengths are very variable (Li et al. 2001). These neurons are the counterparts of eINs in our model. In our simple anatomical model, we therefore added 500 μm ascending axon collaterals to the eINs to carry more excitation rostrally. This changed the distribution of axons [Fig. 4(d, f)] and increased their R-C gradient in the spinal region. Adding 500 μm ascending axons to all eINs in the population model resulted in both a sustained alternating rhythm and a R-C sequence of neuron discharges. Ascending axons are therefore crucial for sustained normal operation of the swimming network.

We then investigated the contribution of ascending eIN axons further by restricting their distribution to approximate those found anatomically. First, the probability of each eIN having an ascending axon was reduced to 0.5 over the whole length of the model. We then further restricted ascending axons, still with a probability of 0.5, to only the more rostral eINs, located down to a caudal limit set at 1,700, 1,200, 1,000 and finally 800 μm from the midbrain. All these models with restricted ascending axons still produced sustained neural activity that alternated on the left and right sides. In addition, the normal R-C sequence of firing became established provided eINs located down to caudal limits of 1,000 μm or more had ascending axons. However, when the caudal limit of eINs with ascending axons was 800 μm or less, the R-C sequence of neuron discharges failed shortly after initiation of swimming. In the models that did show stable R-C firing sequences, the time delays in neuron firing at different positions decreased with the decrease in the number of ascending axons and the consequent shallower R-C gradient in excitation that eINs established. Since the critical feature to obtain activity with rostro-caudal delays was a sufficient number of ascending axons to produce a sufficient rostrocaudal gradient in excitation, in what follows we decided for simplicity to use the model where all eINs had a 500 μm ascending axon.

3.2.3 Sensory initiation of activity

The sensory input that initiates activity can influence the longitudinal position where activity starts. If the sensory excitation went to all neurons, motoneuron activity could start at different longitudinal positions in the network. To get activity to start rostrally and progress caudally it was necessary to connect sensory excitation only to neurons in the rostral 1.5 mm of the model. With this pattern of sensory excitation and a 20 ms delay between the sensory EPSP (SEPSP) on the left and right sides, motor activity started reliably at the head end of the model.

3.2.4 Activity in the population model

During swimming the activity of the individual neurons was now similar to that seen with sharp microelectrodes (Fig. 5). The simulated activity was in the frequency range found in the tadpole (12–25 Hz), was stable, and its frequency increased with the level of synaptic excitation from eINs and decreased with the level of inhibition from iINs. More rostral MNs tended to fire earlier than caudal ones on the same side of the model and the R-C longitudinal delays in firing times with different positions were ∼15 ms/mm. Such delays are much longer than in the tadpole and this issue will be addressed below.

Fig. 5
figure 5

Motoneuron (MN) activities one second after initiation of swimming in the full-length model with no ACh and electrical synapses. (a) Alternating spiking of a left and a right side MN at 1.25 mm. (b) Activity of three left side MNs in different positions [Rostral (R) 0.92 mm, middle (M) 1.49 mm and caudal (C) 2.1 mm] to show rostro-caudal delay indicated by dashed lines marking firing times. (c, d) Three cycles of MN discharge times vs. positions on left side in model (c) when all excitatory interneurons have ascending axons and (d) when only half of the excitatory interneurons with positions more rostral than 1.0 mm have ascending axons. (e) Enlarged view of MN discharge times vs. positions in a single cycle in the spinal region of the model (caudal to position 0.9 mm) with fitted regression line. R-C delays in neuron discharges were calculated by fitting linear regression lines to scatter plots in the spinal region in a cycle by cycle basis and taking the average of slopes of the regression lines in five consecutive cycles

3.3 Analysis of population model function

3.3.1 The distribution of synaptic drive received by motoneurons

What the experiments say

Several previous studies have investigated possible gradients in synaptic drive during swimming. Tunstall and Roberts (1994) made intracellular recordings from putative spinal motoneurons in the tadpole during fictive swimming. They showed that there is a rostral to caudal decrease in the size of the tonic excitatory synaptic input. Regarding inhibition, they found a similar decline in the amplitude of the spike after-hyperpolarization and in the size of mid-cycle inhibition. These findings suggested a R-C gradient in the excitatory and inhibitory synaptic drives received by MNs during swimming. The distribution of excitatory synaptic drive was further investigated by Zhao et. al. (1998) when the non-linear summation of EPSPs was also taken into account. This experimental study showed R-C gradients in glutamatergic excitation but no gradient in cholinergic excitation. A cellular model of the tadpole neurons was used to estimate the synaptic conductances received by MNs on the basis of membrane potential amplitudes. The results suggested that the composition of the synaptic drive to motoneurons changed with longitudinal position. Yoshida et. al. (1998) created a simple spatiotemporal model and used the distribution of the number of axons and the shape of inhibitory synaptic potentials to estimate the spread of inhibition along the cord. They found that the inhibitory synaptic conductance during swimming declined caudally. We have therefore investigated the significance of the longitudinal distributions of these excitatory and inhibitory synaptic drives to MNs using our population model.

What the model shows

We investigated longitudinal gradients in our full-length, 3.5 mm population model with all the connections according to Fig. 3 to see if they matched the synaptic drives found experimentally. First, we simply calculated the distribution of the number of each type of synapse received by MNs. This was done by counting the number of glutamate, ACh, glycine and electrical synaptic connections received by each MN on one side of the model (Fig. 6) These distributions showed clear R-C gradients in the number of glutamate and glycine synaptic contacts but no obvious gradient for ACh and electrical contacts in the spinal region (caudal to 0.9 mm). The presence or lack of R-C gradients in the number of synapses along the cord is related to the distributions of soma and axon lengths of cells, all of which are different for the INs and MNs in our model, as in reality.

Fig. 6
figure 6

Longitudinal distributions of synaptic inputs to motoneurons in the full-length model. The number of glutamate (a), glycine (b), acetylcholine (c) and electrical (d) synapses received by motoneurons at different positions along the R-C body axis on one side of the population model

There is no direct anatomical data on the number of functional synaptic contacts received by spinal neurons in the tadpole. Furthermore, the number of contacts is dependent not only on the distributions of neurons and axons but also on the probability of synapse formation between axons and dendrites. Moreover, physiological features, such as axonal conduction and synaptic delays at synapses, and therefore summation, may alter the distributions of synaptic drives significantly relative to the distributions of synaptic contacts expected on the basis of neuron and axon distributions. To assess gradients, we therefore measured the actual distributions of chemical synaptic drives in the model during swimming activity. We chose 15 left side MNs from rostral, middle and caudal regions of the spinal part of our model. We selected these regions in the model to correspond to those recorded in experiments on the distribution of synaptic drives (Zhao et al. 1998). Since synapses were modelled as conductance changes, the synaptic conductance was dependent on the phase of activity during swimming. We therefore recorded how synaptic conductances changed over time and found the averaged peak conductance during swimming in the selected MNs for each of the three spinal regions.

The distributions of the peak synaptic conductances in MNs during swimming (Fig. 7) showed a clear R-C gradient in glutamate mediated excitation (ANOVA, p < 10−5, slope of linear regression = −13.8 nS/mm, r 2 = 0.79) and in glycine inhibition (ANOVA, p < 3 × 10−5, slope of linear regression = −97.8 nS/mm, r 2 = 0.76), while the ACh synaptic drive showed no clear longitudinal trend (ANOVA, p > 0.05, slope of linear regression = 2.2 nS/mm, r 2 = 0.30).

Fig. 7
figure 7

Distributions of glutamate (Glut), acetylcholine (ACh) and inhibitory (Inh) synaptic drives during swimming in the full-length model. The spinal cord was divided into rostral, middle and caudal regions with coordinates 1.00–1.31, 1.39–1.70 and 1.77–2.34 mm respectively. Five MNs were randomly selected from each region. Maximum synaptic conductances during five cycles 1 s after the initiation of swimming were measured in left side MNs and averaged over each region. Drives were graphed as percentage values taking them as 100% for the region where the MNs received the largest amount of synaptic conductance during swimming. Closed, hatched and open bars are for rostral, middle and caudal regions. Error bars are SDs

3.3.2 Gradient in spatial distribution of excitatory interneurons is necessary to get rostro-caudal delays in motoneuron discharges

Anatomical evidence suggests that both excitatory and inhibitory interneurons are always more densely packed in the rostral than in the caudal spinal cord (Hartenstein 1993). We wanted to check the possible significance of these non-homogeneous neuron distributions found in the spinal cord. To simplify these tests we reduced the length of the model to represent the more rostral region of the spinal cord that has been studied the most physiologically. In our 3.5 mm long model, we chose the region from 1.0 to 2.5 mm which represented the rostral to middle part of the living spinal cord (see Fig. 3). This is the region where there is a clear R-C gradient in excitation to motoneurons [Fig. 6(a)]. During our simulations, the full length model was always built first. We then made cuts at 1.0 and 2.5 mm so the axons connecting neurons in this region to the more rostral and more caudal parts of the full model were cut, just as in physiological experiments with transected spinal cords. In what follows we will describe this model and will refer to it as ‘reduced-length model’ in the figure legends.

In our tests on the shorter, reduced-length, model network we gradually removed the longitudinal gradients in distributions of neurons to reveal if these gradients are responsible for R-C progression of motoneuron discharge. To ensure that effects were due to the distributions of neurons, we always kept the total numbers of neurons the same (57 eINs and 93 iINs per side) when eINs and/or iINs were redistributed. We did not change the distribution or numbers of MNs and their cholinergic and electrical synapses were turned off. We redistributed interneurons in three different ways: we flattened the distribution of either the eINs or iINs without altering the normal R-C distribution of the other type of interneuron, or we redistributed both eINs and iINs evenly. We initiated swimming and checked if: sustained neural activity outlasted the duration of the excitation resulting from the sensory stimulus, firing of neurons alternated between the two sides, and the normal rostral to caudal sequence of neuron discharges was established. We compared these features to the normal case when the R-C gradients in the distributions of these INs were present.

In all three cases, when the distribution of INs was altered, the network activity was persistent and the firing of left and right side INs and MNs alternated. However, the size and the sign (R-C vs. C-R) of the delays in firing of neurons with different axial positions were affected (Fig. 8). In the control case with the unaltered distributions of INs, the R-C delay was about 14–15 ms/mm, higher than the experimental values [Fig. 8(a), but see below]. When the distribution of iINs was flattened but the normal gradient in distribution of eINs was kept [Fig. 8(b)] the R-C firing sequence remained and the size of delays did not change significantly. On the other hand, when the eIN R-C gradient was flattened but the R-C gradient kept for the iINs [Fig. 8(c)], the sequence of motoneuron discharges was reversed and the caudal neurons fired before more rostral ones. A similar reversal of firing sequence was obvious in most cycles when both the eINs and iINs were given flat distributions [Fig. 8(d)]. However, in this case the sign of the delay was changeable, sometimes without significant R-C trends in firings during the cycles, and linear regression lines had the lowest correlations. Our simulations suggest that it is the R-C gradient in spatial distribution of excitatory interneurons that make rostral motoneurons fire prior to more caudal ones, as seen in normal forward swimming.

Fig. 8
figure 8

Effect of changed distributions of interneurons on longitudinal coordination in the reduced-length model with no ACh and electrical synapses. (a–d) Show discharge times vs. position for left side motoneurons in models with different distributions of interneurons along the spinal cord. Times of spikes are shown in three consecutive swimming cycles 1 s after initiation of network activity. Icons show if distributions of excitatory (eIN) and inhibitory (iIN) neurons along the R-C body axis decrease caudally or remain constant. Possible R-C trends in neuronal discharges were tested by the slope of linear regression lines. (a) normal case with R-C gradients in distributions of eINs and iINs (R-C delay, 14.2 ± 1.3 ms/mm, p < 10−22, r 2 = 0.73–0.80). (b) Gradient in eINs only (R-C delay, 16.5 ± 0.5 ms/mm, p < 10−24, r 2 = 0.77–0.81). (c) Gradient in iINs only, (R-C delay, −4.6 ± 0.1 ms/mm, p < 10−4, r 2 = 0.21–0.23). (d) Flat distribution for both eINs and iINs (R-C delay, −2.6 ± 3.0 ms/mm, p = 10−8–0.42, r 2 = 0.01–0.37)

3.3.3 Motoneuron to interneuron ACh synapses decrease the rostro-caudal delays in motoneuron discharges and synchronise motoneurons locally

Perrins and Roberts (1995a, b, c) made simultaneous intracellular recordings from pairs of MNs and revealed that they are interconnected by chemical and electrical synapses. Moreover, they proposed that MNs make feedback cholinergic synaptic contacts onto premotor interneurons of the spinal CPG (Perrins and Roberts 1995c, Fig. 3, but see Section 4). This type of connection would make MNs members of the CPG, rather than simple output elements driven by the premotor INs and exciting the muscles. On the other hand, rhythm generation is still possible in the tadpole even when ACh mediated synaptic transmission, is blocked with the nicotinic ACh receptor antagonist DHβE (Zhao et al. 1998). This antagonist caused a small decrease in excitation of rhythmic neurons, and reduced the reliability of spike production. These results suggest that central cholinergic synapses may modulate rather than play an essential role in rhythm generation.

We investigated the possible effects of the MN-to-IN feedback synapses in the shortened model (1.0 to 2.5 mm). With the normal distributions of INs and with cholinergic connections between MNs, we looked at the consequence of adding the motoneuron-to-interneuron feedback ACh synaptic connections on the R-C delays. The timing of motoneuron discharges were sampled in rostral, middle and caudal segments on one side of the cord, with and without the feedback MN connections added (Fig. 9). With the motoneuron-to-interneuron feedback connections turned off the R-C delay was much too high (29.2 ± 6.3 ms/mm) relative to the physiological values (2–5 ms/mm; Tunstall and Roberts 1991). When feedback connections were added, the R-C delay decreased considerably to 7.5 ± 0.7 ms/mm (t test, p < 0.003) nearing the physiological range of delays [Fig. 9(a)].

Fig. 9
figure 9

Rostro-caudal delays and ventral root burst durations with and without motoneuron-to-interneuron feedback connections in the reduced-length model with no electrical synapses. (a) Discharge times vs. positions of MNs in three spinal segments were plotted in a cycle, 600 ms after initiation of swimming. Slopes of linear regression lines revealed that R-C delays dropped significantly (t test, p < 0.003) to near physiological levels when motoneuron feedback connections were added. Data points and equations of linear regression lines are: closed circles, y = 0.0067x + 0.7194 (r 2 = 0.48) and open circles, y = 0.0367x + 0.6191 (r 2 = 0.87) with and without MN-to-IN connections respectively. (b) Ventral root burst durations in three spinal segments measured as the maximum difference between times of spikes in MNs within the same segments. Closed bars are for burst durations with the feedback connections and open bars without. Error bars are SDs

The scatterograms and r 2 values also suggested that feedback ACh synapses not only decreased the R-C delays but also synchronised the motoneuron discharges within the same spinal segment. To investigate this further, we approximated the duration of ventral root discharges as the maximum difference in firing times of motoneurons located within the same segment. Comparison of ventral root discharge durations [Fig. 9(b)] showed a decrease from 10.9 ± 1.6, 25.4 ± 2.5, 15.2 ± 4.9 ms to 5.6 ± 0.0, 7.5 ± 0.5, 13.6 ± 1.2 ms in the rostral, middle and caudal segments respectively if feedback ACh connections were added. These durations are in the physiological range (5–10 ms, Dale 2003) except in the caudal segment (but see below).

3.3.4 Electrical synapses synchronise motoneuron discharges locally

Perrins and Roberts (1995a, b) have shown that neighbouring spinal MNs are interconnected by electrical synapses in the tadpole. However, in the absence of drugs which are specific enough to block electrical synapses (see: Rozental et al. 2001; Rouach et al. 2003; Leznik and Llinas 2005; Srinivas and Spray 2003; Wang et al. 2006), it is difficult to assess the significance of electrical coupling in the spinal CPG for swimming in the whole tadpole. Therefore we conducted simulations to compare the patterns of activities with and without motoneuron electrical synapses and by varying the coupling strengths of these connections.

The junctional to resting neuron conductance ratios have been estimated to be between 10% and 23% in spinal MNs of the tadpole during swimming and during stimulation of sensory pathways (Zhao et al. 1998). Taking the resting neuron conductance as 1 nS for a non-shunted MN, leads to an estimate of 0.1–0.2 nS (5,000–10,000 MΩ) for the junctional conductances. Perrins and Roberts (1995a) supplied experimental values between 0.1 and 0.2 for the coupling ratios between MNs in the tadpole spinal cord. Using the formula for the relationship between the coupling ratio and the leakage conductance of the neuron derived by Bennett (1966) the unitary synaptic conductance for the electrical synapses yields values of 0.11–0.25 nS (4,000–9,000 MΩ). Galarreta and Hestrin (2002) have also estimated a similar 5,000 MΩ mean coupling resistance between pairs of interneurons in the mouse neocortex.

Therefore initially, electrical couplings of 5,000 MΩ were added between all MNs within each 150 μm long segment. MNs located in different segments were not interconnected by electrical synapses but ACh synapses from rostral to more caudal MNs and the MN-to-interneuron feedback synapses were modelled. When electrical synapses were added, the duration of ventral root discharges decreased from 5.6–13.6 to 3.5–6.3 ms but the R-C delays in MN firing times increased slightly (from 7.5 ± 0.7 to 9.3 ± 1.3 ms/mm, t test, p = 0.045) [Fig. 10(a)]. This brings the durations of ventral root bursts in the simulated CPG below the upper limit of the physiological range of 5–10 ms (Dale 2003). We then varied the resistance of electrical couplings between infinity (no coupling) and 100 MΩ. Generally, the duration of ventral root discharges decreased with increasing coupling strength (decreasing coupling resistance) indicating better intrasegmental synchronisation of motoneuron activity. To measure the dispersion of ventral root burst durations along the cord, the ratio of the maximum and minimum durations among ventral roots were calculated. We found that with strengthening electrical coupling, the differences among burst durations in ventral roots of various segments was increased [Fig. 10(b)].

Fig. 10
figure 10

The effects of electrical coupling between motoneurons on intrasegmental synchronization of motoneuron discharges and on dispersion of ventral root burst durations along the spinal cord in the reduced-length model. (a) Ventral root burst durations in a rostral, middle and caudal spinal segment with different sizes of electrical couplings and with no coupling between MNs. Positions of segments were the same as in Fig. 7. Closed, hatched, dotted and open bars are for no electrical coupling, 5000, 500, and 100 MΩ coupling resistances respectively. (b) Relative burst durations as a function of strength of electrical coupling between motoneurons. Relative burst durations were calculated as the ratios of the maximum and minimum burst durations of the ventral roots. Error bars are SDs

3.4 How critically are the network properties dependent on the properties of neurons?

To investigate the generality of our findings, we created a new model neuron with features that correspond to spinal rhythmic neurons recorded by patch electrodes in the tadpole (Li et al. 2004a). The new model neurons had a less negative resting membrane potential (−55 mV), spikes that overshoot zero potential, and fired a single spike in response to current injection like spinal excitatory interneurons (Li et al. 2006; model neuron type 2, see Appendix for full list of features). We replaced the neurons in our network with the new one and tested if network properties changed as a result of this change in neuron properties. By just increasing the unitary glutamate and inhibitory synaptic conductances to 0.5 and 10 nS, the model produced stable swimming in the reduced-length spinal and in the full-length model with and without ACh and electrical synapses. The frequency range and the R-C delays of MN discharges were smaller (14–16 Hz and 12–15 ms/mm). Having checked the new network model, we systematically repeated our previous simulations and compared the results.

In simulations with the first type of model neuron, which mimicked CPG neurons recorded by sharp microelectrode, the gradient in distribution of eINs proved to be necessary for the network to generate the R-C sequence in MN firing. We found the same with type 2 neurons. The sign of the delay (rostro-caudal vs. caudo-rostral) was determined by the presence or absence of a gradient in the distribution of eINs along the spinal cord. We then looked for the effect of MN-to-IN feedback connections. When these connections were added, both the R-C delays and durations of ventral root discharges tended to decrease, as in the previous simulations. With type 2 neurons the durations of ventral root bursts were shorter, both with and without these feedback connections, especially in the rostral and middle spinal segments. With the feedback connections added, burst durations were within the physiological range (below 10 ms) having the shortest values in the most rostral region of the cord. When 5,000 MΩ electrical connections between MNs within segments were also added, the lengths of bursts decreased (from 3.1 ± 0.2 to 2.8 ± 0.1 ms: rostral, 6.2 ± 0.9 to 5.8 ± 1.2 ms: middle, and 8.0 ± 1.2 to 7.3 ± 1.3 ms: caudal ventral roots). As electrical coupling was strengthened, burst durations decreased by a few percent in all parts of the cord. However, when electrical coupling was strong, 100 MΩ, differences between the ventral root burst durations along the cord increased.

Generally, the replacement of the model neuron did not change any of our previous conclusions; it only resulted in minor shifts in quantitative measures of the patterns of activity generated by the swimming network suggesting a relative independence of our findings on the detailed cellular properties.

4 Discussion

To study the longitudinal and local coordination in the spinal CPG for swimming of the tadpole, we created a computational model. The model had a realistic R-C length and used anatomically realistic densities of neurons and axons lengths, which were, for example, dependent on the position of neuron cell bodies. Rules of network connections were based on detailed knowledge about the CPG from anatomy, pharmacology and physiology.

4.1 Significance of distributions of interneurons and axons along the spinal cord

We have made use of this population model to study the significance of the distribution of premotor interneurons along the spinal cord. We have shown that a R-C gradient in the density of excitatory interneurons is necessary for the normal R-C sequence of motoneuron firing times seen during swimming. This produces a R-C wave in contraction of the swimming muscles, and a mechanical wave of bending along the trunk which propels the animal forward in water. We also showed that the normal fall in the number of interneurons along the length of the spinal cord, together with the connection rules adopted, resulted in distributions of Glutamate and Glycine synapses with R-C gradients along the cord. We suggest that this R-C gradient depends on an underlying gradient in axon and therefore synapse numbers and illustrates the way in which axon numbers depend on neuron distribution and axon length and projection direction (Fig. 4). Since there is not always a simple relationship between the distribution of the number of synapses made by interneurons and the distribution of synaptic drive to motoneurons during swimming (Li et al. 2001), we checked the R-C distributions of excitatory and inhibitory synaptic drives (summed synaptic conductances) in our model against the ones derived earlier (Zhao et al. 1998; Tunstall and Roberts 1994; Yoshida et al. 1998). We concluded that these synaptic drives were very similar to those calculated previously. Therefore, our model suggests that the experimentally observed R-C gradients in the distributions of excitatory interneurons, in the excitatory synaptic drive to motoneurons and in the firing times of motoneurons, are closely interrelated with the production of normal forward swimming.

4.2 Functional roles of the presumed excitatory motoneuron-to-interneuron feedback synapses

The feedback connections from motoneurons to interneurons were proposed by (Perrins and Roberts 1995c) on an experimental basis. They found that interneurons receive rhythmic cholinergic excitation during swimming from within the spinal cord that could be blocked by nicotinic receptor antagonist. Since at the time of these experiments the only known type of rhythmically active cholinergic neuron of the spinal CPG for swimming in the Xenopus embryo was the motoneuron, they proposed them as the most likely source of this excitation. More recently, the discovery of co-release of Glutamate and ACh at synapses made by the descending interneurons (Li et al. 2004c) put this interpretation into doubt but did not rule out the option of motoneuron-to-interneuron feedback connections. Therefore, we studied the possible roles, if any, of such connections during swimming.

Our population model suggests that motoneuron-to-interneuron feedback synapses decrease the intersegmental lag in firing times, reducing the R-C delays. This is likely because MN axons are more numerous caudally since they have only descending axons, and in this way their synapses can excite more caudal neurons to fire earlier than without these connections. Apart from reducing the R-C delays, we also found that these feedback synapses can enhance the local synchronizing effects of electrical synapses between motoneurons. Thus we could attribute two possible roles of these feedback connections and support the original proposal of Perrins and Roberts (1995c) that would make MNs parts of the spinal CPG for swimming.

4.3 Electrical synapses between motoneurons

Electrical coupling among pattern-generating elements is common (Kiehn and Tresch 2002; Connors and Long 2004). It has been long appreciated that electrical coupling tends to synchronize the activity of neurons. However, there are examples where coupled neurons do not fire synchronously and may even fire out of phase (Sharp et. al. 1992). Sherman and Rinzel (1992) have shown theoretically that out of phase activity can be produced if coupling strength is low (see also Marder and Calabrese 1996). Cellular and membrane properties may also shape the effect of electrical coupling.

Our simulations confirmed that the intrasegmental synchrony in firing of motoneurons was improved by increasing the coupling conductance between motoneurons. However, we found that a compromise is needed between the dispersion of ventral root burst durations along the spinal cord and the degree of intrasegmental synchronization of motoneuron discharges: synchrony is improved but differences among burst durations are increasing with tighter coupling.

4.4 How critically are our conclusions dependent on the cellular properties of neurons?

For our simulations we used model neurons with the same membrane properties and all fired just once to current injections. However, we built the population model of the swimming CPG with two different model neurons, without modifying the distributions of neurons and axon lengths or altering the principles of synaptic connectivity in the default population model. Although changing the model neurons in simulations caused some changes in the size of R-C delays per unit length along the cord and in the duration of ventral root discharges, these changes in the absolute values did not alter the general conclusions on the significance of motoneuron-to-interneuron feedback connections and on the effects of electrical couplings among motoneurons. The R-C distributions of synaptic drives to motoneurons were also virtually unaffected by the changed membrane properties of the model neurons. All these features suggest that our major observations on the significance of interneuron distributions, on the proposed feedback synapses to interneurons, and on the interrelationship between the distributions of synaptic drives to motoneurons and the normal R-C sequence of neuron discharges are all rather insensitive to detailed membrane properties. Of course, we do not want to underestimate the importance of neuron specific properties in the spinal CPG (Sautois et al. 2007; Li et al. 2007a). On the other hand, in a recent simulation study Sautois et al. (2007) compared the swimming CPG of the tadpole when all neurons had the same membrane properties and when neurons of the network were given cell specific features. They found that the CPG with identical neurons with the property of producing only single spikes to current injections produced similarly stable swimming pattern with the same basic features as the CPG with different more realistic cellular properties.

This relative independence from membrane properties may be interpreted as the idea that behaviorally crucial features, vital for survival, like normal forward swimming, must be based on very simple developmental rules, like linear gradients in the density of interneurons along the spinal cord (Li et al. 2007b).

4.5 Comparisons with lamprey swim models

Many of the issues addressed in this study on young Xenopus tadpoles have also been addressed in extensive modelling studies on another swimming animal: the lamprey (reviewed in: Grillner et al. 2007). In both cases a CPG coordinates alternating contractions of swimming muscles which pass from head to tail and attempts have been made to construct models of the neuronal circuitry responsible. Some of the lamprey modelling has been “top down”, applying coupled oscillator theories and tested predictions using experiments (Cohen et al. 1992). More often the modelling strategy has been “bottom up”: starting with experimental data on membranes, neurons and synapses and working up towards the construction of whole CPG neuronal networks. However, vertebrate CNS circuits are sufficiently complex to make this approach less successful than it has been in simpler systems like the circuits controlling the longitudinal spread of the leech heartbeat (Jezzini et al. 2004).

In the lamprey, “bottom up” modelling is based on good evidence on the properties of spinal motoneurons that allow them to generate bursts of activity when continuously excited by bath application of NMDA (Huss et al. 2007, 2008). Mutual excitatory connections between spinal premotor excitatory interneurons form the basis of recent lamprey studies (e.g. Kozlov et al. 2007). The first proposal on the importance of such connections came in the tadpole and was based on indirect evidence (Dale and Roberts 1985). They were then incorporated into early models of the tadpole spinal swimming circuit (Roberts and Tunstall 1990) but direct evidence for the existence of these synapses in the tadpole has only been obtained recently (Li et al. 2006). This direct evidence for mutual excitatory connections between premotor interneurons provides one of the bases for our modelling study. We have argued that in the tadpole such connections are essential if swimming activity is to self-sustain after a brief sensory stimulus (Li et al. 2006). Despite their appearance in diagrams and importance in recent models (Lansner et al. 1997; Grillner et al. 2007; Kozlov et al. 2007; Várkonyi et al. 2008), there is still no evidence for equivalent mutual excitatory connections in the lamprey (Grillner et al. 2007).

In lamprey spinal cord models, activity has generally been evoked, as in experiments, by bath application of excitants like NMDA and is usually not self-sustaining. This is a significant difference to the tadpole where long-lasting swimming is evoked by a brief skin stimulus in preparations that have to include at least some hindbrain (Li et al. 2006). The term CPG is applied to both types of network but in one the CPG may simply organize alternating activity when provided with continuous excitation whereas in the other (tadpole) the CPG has a built-in mechanism to sustain its own activity without continuous external excitation.

In both lamprey and tadpole there is a head-to-tail spread of motor activity during swimming. We have shown here how the underlying head-to-tail gradient in excitation that is seen during swimming in the tadpole (Tunstall and Roberts 1994; Zhao et al. 1998) may result from the higher rostral density of excitatory interneurons when coupled with their rather short axonal projection patterns (Li et al. 2001, 2006). This organisation leads to a gradient in axon numbers and our modelling shows that this may be sufficient to account for a longitudinal gradient in synapses and excitation. In multisegment lamprey models, neurons have usually been evenly distributed, as in the model by Wadden et al. (1997), who used a model without segmental boundaries. If the projection distances and synaptic strengths were suitably tuned, this model produced alternating motor activity which spread from head-to-tail based on the stronger caudal reciprocal inhibition but the evidence for this in the animal is unclear. Later models gave excitatory interneurons longer descending than ascending axonal projections, based on experimental evidence (Dale 1986), and showed that rhythmic bursting activity could become established and progress from head to tail, even within a single half of the spinal cord (Kotaleski et al. 1999; Kozlov et al. 2007).

Besides the role of a R-C gradient in the number of excitatory interneurons along the cord, the present investigation has also suggested that the sensory excitation that starts swimming needs to be stronger rostrally if the activity is to start at the head end. At present we have no physiological evidence for such a gradient, but it is expected because both the primary sensory RB neurons and the sensory relay interneurons that they excite project rostrally where their axon numbers will build up and lead to stronger excitation (Clarke et al. 1984; Li et al. 2001).

We hope that this discussion and comparison of some selected and relevant issues will have pointed to some of the difficulties in attempts to construct evidence based models of complex CNS circuits and how these difficulties increase as one moves from simpler (leech heart) to more complex (lamprey) circuits. In all studies it remains essential to distinguish where the evidence is not yet available.

4.6 Final remarks

Without a continuous population model with R-C length mimicking real anatomy, it would not have been possible to dissect the significance of the spatial distributions of neurons, the presumed cholinergic motoneuron-to-interneuron feedback or the controversial effects of electrical synapses. Although many of the predictions by the model cannot be yet tested experimentally, our model has increased our understanding on the inter- and intrasegmental coordination in the tadpole CPG during swimming. On the other hand, in the absence of relevant experimental information we had to make several assumptions. While we have acceptable information on distributions of neurons and their projection distances, in reality we need much better data, for example based on the transcription factors expressed by each different type of spinal neuron. We also assumed that neurons and synapses of given kinds have exactly the same properties, and that synaptic contacts between neurons and passing axons have the same probabilities of formation along the cord. We know that this is not correct and that there is variability. We are beginning to get information on synapse contact probabilities and these can now be used in future modelling (Li et al. 2007a).