Abstract
How do the pioneer networks in the axial core of the vertebrate nervous system first develop? Fundamental to understanding any full-scale neuronal network is knowledge of the constituent neurons, their properties, synaptic interconnections, and normal activity. Our novel strategy uses basic developmental rules to generate model networks that retain individual neuron and synapse resolution and are capable of reproducing correct, whole animal responses. We apply our developmental strategy to young Xenopus tadpoles, whose brainstem and spinal cord share a core vertebrate plan, but at a tractable complexity. Following detailed anatomical and physiological measurements to complete a descriptive library of each type of spinal neuron, we build models of their axon growth controlled by simple chemical gradients and physical barriers. By adding dendrites and allowing probabilistic formation of synaptic connections, we reconstruct network connectivity among up to 2000 neurons. When the resulting “network” is populated by model neurons and synapses, with properties based on physiology, it can respond to sensory stimulation by mimicking tadpole swimming behavior. This functioning model represents the most complete reconstruction of a vertebrate neuronal network that can reproduce the complex, rhythmic behavior of a whole animal. The findings validate our novel developmental strategy for generating realistic networks with individual neuron- and synapse-level resolution. We use it to demonstrate how early functional neuronal connectivity and behavior may in life result from simple developmental “rules,” which lay out a scaffold for the vertebrate CNS without specific neuron-to-neuron recognition.
Introduction
In vertebrates, molecular mechanisms, like gradients of morphogens, organize dorsoventral (DV) rows of neurons along the developing nervous system (Helms and Johnson, 2003; Goulding and Pfaff, 2005; Lewis, 2006). The neurons then grow axons, and knowledge of the mechanisms allowing axons to navigate is constantly increasing (Dickson, 2002; Chilton, 2006; Zou and Lyuksyutova, 2007). However, it is still unclear how neurons build networks by making synaptic connections (Sperry, 1963; Zipursky and Sanes, 2010). Our aim is to “grow” realistic model networks to test whether early network assembly could be controlled by a basic set of rules. Critically, to validate model networks we need to show they have the same responses to input as living networks. This raises a big problem for exploring networks in higher brain regions because precise inputs and outputs are not known. In newly hatched larval fish and amphibians, pioneer networks allow them to respond predictably to touch by swimming (McLean and Fetcho, 2009; Roberts et al., 2010). The precise input and output are therefore defined and can be used to validate model networks.
We use hatchling Xenopus tadpoles to explore the first formation of a working network (Fig. 1A–C; Roberts et al., 2010). Paired whole-cell recordings have provided detailed evidence on the different spinal neurons and synaptic interactions in the network generating swimming in response to touch stimuli. This evidence provides the foundation for using a developmental model to generate a full-scale functioning locomotor network. The recordings and previous anatomical network modeling suggested a lack of specificity in synaptic connections (Li et al., 2007). Connection probability appeared to reflect the dorsoventral distributions of axons and dendrites (Borisyuk et al., 2008, 2011): dorsal sensory axons contact dorsal sensory pathway neuron dendrites, but not ventral motoneuron dendrites. Our hypothesis is that early axon growth is controlled by simple responses to basic features like physical barriers and chemical gradient cues (Sperry, 1963; Shimozono et al., 2013); axons synapse probabilistically with any dendrites they contact without specific neuron recognition; this is sufficient for networks to form that are able to function without further refinement.
If the factors controlling network assembly in the early vertebrate brain are simple, a developmental model using basic rules should be sufficient to generate the full-scale synaptic connectivity of the functional network controlling swimming in Xenopus tadpoles. We therefore use existing (for review, see Roberts et al., 2010) and new whole-cell recordings to define neuron activity and obtain anatomical data on the seven neuron types involved. We then build an axon growth model to reconstruct the swimming network of ∼1400 neurons in the tadpole hindbrain and rostral spinal cord. Crucially, we know exactly what sensory activity induces tadpole swimming as well as the resulting motor output. This means we can test the operation of networks generated by our developmental model by asking if they produce swimming activity in response to sensory stimuli. We conclude that simple mechanisms may lay out the first functioning neuronal networks in the vertebrate nervous system.
Materials and Methods
Animals, behavior, physiology and anatomy
Procedures for obtaining hatchling Xenopus laevis (Daudin) tadpoles of either sex from a captive breeding colony comply with UK Home Office regulations. All experiments and analysis were approved by the local ethical committee and were performed on tadpoles at developmental stage 37/38 (Nieuwkoop and Faber, 1956) at 18–22°C.
Videos of 40 tadpoles were recorded at 300 fps with a Casio Exilim ExF1 camera. Methods for extracellular skin stimulation, motor nerve recording, whole-cell electrophysiology, and dye filling with neurobiotin have been described recently (Buhl et al., 2012). Stable whole-cell recordings were made in 74 animals: 77 descending interneurons (dINs; including six pairs simultaneously on opposite sides), 22 dorsolateral commissural interneurons (dlcs; including 11 simultaneously with a dIN), and 3 dorsolateral ascending interneurons (dlas). Following experiments, neuron anatomy was revealed using DAB as chromogen and observed using a 100× oil-immersion lens (for review, see Li et al., 2001). The recorded dlc, dla, and dIN neurons were identified by comparison to published data on their responses to skin stimulation, activity during swimming, and anatomy after dye-filling (Li et al., 2003, 2004a, 2006). Chemicals were from Sigma.
To measure the three-dimensional morphology of neurons, neurobiotin-filled neurons in CNS whole mounts (Fig. 1F) were viewed at 20× or 40× using a modified Nikon Optiphot microscope with cool LED illumination (Cairn) using a DeltaPix DP 200 camera. Specimens were positioned in three dimensions using a computer-controlled Scientifica Patchstar manipulator in place of the microscope stage. Using the Patchstar Gridstore software, neuron structures were lined up against cross-hairs and pairs of 3D coordinates were recorded (resolution 1 μm) defining the diameter of neuron segments, perpendicular to their long axis. The cross-hairs were also used to define the positions of the vertical hindbrain–midbrain border and the horizontal CNS ventral midline. These “axes” and their intersection provided datum lines and a CNS framework in which to map neuron positions and features. The type of structures and their 3D coordinates were saved to a text file by Gridstore software. The saved 3D x, y, z coordinate pairs were then corrected for z-axis distortion by light refraction in glass (1.35×), shrinkage (1.28×), and straightening of curvature of the ventral midline “axis.” The x, y, z coordinate pairs were converted to .swc format with a single coordinate at the midpoint of the pair plus a radius. Transformations were performed using custom Python and MatLab tools (available on request). Swc files list x, y, z midpoints and their radii along with connectivity and type, and allow 3D visualization of neurons based on a series of connected cylindrical segments using Neuromantic (Myatt et al., 2012). This also allowed the export of images (Fig. 2D–F). For use in the growth model (see below), swc coordinate data were converted to 2D as if viewed laterally, and the dorsoventral projection patterns of dendrites and axons were defined for each neuron type relative to the ventral midline axis (Fig. 2F,G). The new swc neuron mapping data (79 neurons; including 8 dIN fills from recordings made during this study) were combined with new measurements made from existing drawing tube tracings of individual neurobiotin or HRP-filled neurons (205 neurons). These new measures extended and improved the accuracy and detail of our anatomical library. We measured soma and axon branching positions, axon initial outgrowth angle and tortuosities, dorsoventral dendrite extents, and complete primary and secondary axon trajectories and lengths. Calculations and descriptive statistics were performed in Microsoft Excel and Minitab. All means are given with their SD.
To define CNS anatomy and the axon growth environment, tadpoles were anesthetized in 0.1% MS222, fixed in 2.5% glutaraldehyde in 0.1 m cacodylate buffer, pH 7.4, at 4°C for 90 min, and transferred to buffer to make direct measurements of whole-mount CNS dimensions before dehydration (n = 4). Three specimens were postfixed in 2% osmium tetroxide, dehydrated, embedded in Epon 812, sectioned at 2 μm, and stained with 1% methylene blue. Photographs of sections (Fig. 2A) were used to measure CNS features using NIH ImageJ software.
Axon growth and synapse formation model
On the basis of our measurements, we open the CNS along its dorsal midline like a book to produce a two-dimensional growth environment representing the outer part of the tadpole CNS where the axons grow (Borisyuk et al., 2011). The axon growth model is described by the following system of three difference equations: where (xn, yn) are the variables defining the coordinates of the axon tip (growth cone) in a two-dimensional field (with a longitudinal x-axis and the y-axis representing the DV axis on both sides of the CNS) and where the variable θn is the growth angle. Parameters are as follows: Δ is the growth step (1 μm); ζn is a random variable in the range [−α, α]; gR, gV, gD are sensitivities to rostral, ventral, and dorsal guidance cues. Guidance cues hR, hV, hD are given by their “diffusive” rostral, ventral, and dorsal gradients: where βR, βV, βD are slopes. All parameters are universal for all neuron types except for the three gradient sensitivities and random angle fluctuations: (gR, gV, gD, α). Values of these four parameters were selected for each axon type using a stochastic optimization algorithm (“patternsearch” routine from the Global Optimization Toolbox of MATLAB) to minimize a cost function. This cost function was given by the weighted sum of two components: (1) the sum of squared differences of the averaged tortuosities (arc-chord ratios) calculated for the generated and measured axons, and (2) the dissimilarity (“distance”) between the dorsoventral distributions of measured and modeled axon projections (histograms). The distance between histograms of measured and modeled axons is the normalized sum of squared differences between counts in corresponding bins of the histogram, and this quadratic form is identical to that used in the χ2 statistic for testing the homogeneity of two distributions.
The optimization procedure results in parameter estimates that ensure the value of the cost function is suitably small: in the majority of cases, the value after optimization is in the range (2–16). To judge the quality of optimization, we applied two sample t tests to show that differences between mean tortuosities for real and modeled axons were not significant (p > 0.05 in all cases). Although the standard statistical tests are not applicable for estimating the similarity between modeled and experimental histograms of dorsoventral axon distributions (data are not independent, being close successive points on the same axon), in the majority of cases, the values of the χ2-based component of the cost-function lay below the appropriate critical value for the χ2 statistic (16.92 for 9 df at a significance level of p = 0.05).
Modeled dendrites were specified by their dorsal and ventral extents, by adding two-dimensional Gaussian noise (equal SDs of 15 and correlation of 0.8) to pairs of measured values. Dorsoventral projections of modeled dendrites were binned (10 μm bins) to produce a histogram. These binned distributions matched those of measured dendrites for each neuron type (χ2 test: p > 0.05). Methods have been described previously for distributing populations of somata for each different neuron type along the CNS, and for allowing the stochastic formation of en-passant axon-to-dendrite synaptic connections to generate a network (Li et al., 2007; Borisyuk et al., 2011). The axon growth and synapse formation model was implemented using custom-written MatLab code.
Conductance-based neuronal network model
To assess each network, we mapped it onto a physiological, conductance-based, neuronal network model to test whether the connections generated by the growth model could produce swimming-like motoneuron activity in response to brief “sensory” stimuli.
Active channels in model dIN neurons.
Hindbrain and spinal cord dINs are central to the Xenopus tadpole swim Central Pattern Generator (Li et al., 2010; Li, 2011). To construct model dINs with appropriate physiological firing properties, we used NEURON (controlled by a Python script; Carnevale and Hines, 2006) to build networks with 30 multicompartment dINs electrically coupled by gap junctions on their axons. Membrane potential (V) evolves according to Equation 3. Channel kinetics were based on voltage-clamp data for the currents (Dale, 1995; Winlove and Roberts, 2012). The sodium (iNa) and potassium currents (iKf, iKs) were modeled using a Hodgkin–Huxley-type formulation and the calcium current (iCa) was modeled using the Goldman–Hodgkin–Katz formulation. The ilk current is due to passive membrane leak conductance and iext is an externally applied current that is only used for testing the response of neurons to current injections.
Each channel is gated by one or more gating variables. Equation 4 describes the dynamics of an arbitrary gating variable X, where the functions X∞ (V) and τ∞ (V) give its steady-state value and characteristic time constant, respectively (Eq. 5). These functions are dependent on functions αX (V) and βV (V), which control the opening and closing of gates and are written in the form of Equation 6. The values for parameters A, B, C, D, and E are given in Table 1. [Note that the β-rate equation for calcium takes one of two sets of parameters depending on the membrane potential (Dale, 1995).] Parameter sweeps over the channel densities (of lk, Na, Kf, Ks, Ca) were used to match the responses of the electrically coupled, multicompartmental model dINs to those of experimentally recorded dINs. The responses matched were as follows: single action potential to depolarizing current; rebound firing after hyperpolarization during depolarization; rhythmic firing during NMDAR activation (Soffe et al., 2009; Li et al., 2010). Since variability is introduced into model parameters (see below), 10 neurons were evaluated at each point.
Physiological simulation of the full network.
For efficient computational simulations of the entire network with ∼1400 neurons, the multicompartmental dIN model was simplified to a single compartment. The forms of the equations remained the same and values for density parameters were converted into whole-cell values by assuming a surface area of 1000 μm2. (This is a round number lying within the range of summed soma and dendrite surface areas measured for filled neurons: 754–1594 μm2; n = 5.) The network model included electrical coupling with a conductance of 0.2 nS between pairs of dINs within 100 μm of each other. This gave coupling coefficients of 5–10% that decreased with distance as in previous experiments (Li et al., 2009). The additional conductance introduced through gap junctions was compensated for by reducing the leak conductance to 1.4 nS so dIN input resistance was in line with experimentally measured values (Sautois et al., 2007). An electrically coupled population of single compartmental model dINs exhibited the same characteristic responses to current injection and simulated NMDA perfusion as the multicompartmental model dINs.
The other, non-dIN, swim neurons in the tadpole do not have the specialized firing properties seen in the dINs. Their “typical” repetitively firing properties were modeled using the equations and parameters of the model motoneuron described by Sautois et al. (2007). For both the dINs and non-dINs, a capacitance of 10 pF was used, corresponding to a density of 1.0 μF/cm2 on a neuron with a surface area of 1000 μm2. The membrane properties for dIN and non-dIN neurons are given in Table 2.
The synapses in the model were implemented as described previously (Sautois et al., 2007). Glutamatergic synapses were modeled as separate conductances activated by AMPAR or NMDAR (with voltage dependency), but only dIN-to-dIN synapses activated both these conductances. Inhibitory glycinergic synapses were given a reversal potential of −75 mV (Dale, 1995). Following a presynaptic spike, the synaptic opening and closing variables were increased by a step value of 1.25 for glutamatergic synapses and 3.0 for glycinergic (changed from Sautois et al., 2007 to better match physiology). The model also includes two components of delay: a fixed synaptic delay (1 ms) and an axonal conduction delay that depends upon the longitudinal distance between the presynaptic and postsynaptic neurons (0.0035 ms/μm).
Synaptic conductances had values based on paired recordings (Sautois et al., 2007) with the following exceptions: sensory RB neuron AMPAR excitation was made stronger (0.6–8.0 nS) onto sensory pathway neurons to compensate for their reduced input resistance in our simplified model; the AMPA conductance of dIN synapses onto inhibitory ascending interneurons (aINs) was reduced (0.6–0.1 nS) so aIN firing during swimming was restricted to a short period following sensory stimulation in agreement with experimental data (Li et al., 2004b); the NMDA conductance of dIN-to-dIN synapses was reduced (0.3–0.15 nS) so the level of steady depolarization in dINs during swimming matched experimental measurements more closely (Li et al., 2006).
Our neuron network model was implemented with custom-written C code. All numerical simulations were performed using the Runge–Kutta–Felberg method from the GNU Scientific Library (version 1.15) with an adaptive step size (absolute and relative tolerance 10−5). The maximum step size was 0.1 ms and spikes were detected after every step by membrane potential zero crossing. Gaussian noise with a SD equal to a percentage of the mean value was applied to the capacitance and leak/channel conductances of all neurons (at 2% of the mean), as well as the conductances of every individual synapse (at 5% of the mean). This noise represents the nonhomogeneity of soma sizes and synapse strengths, respectively. Higher noise levels were tested and found to make no major qualitative difference to the results of simulations.
The code for all our anatomical and physiological modeling is available on request from R.Borisyuk{at}plymouth.ac.uk
Results
The foundations for this study lie in extensive previous work, first, on the neurons controlling tadpole swimming (Roberts et al., 2010; Li, 2011) and, second, on modeling of the operation of the swimming network (Sautois et al., 2007) and its development (Li et al., 2007; Borisyuk et al., 2008, 2011). Our previous axon growth model was not based closely on biology but provided a mathematically simple way to generate the longitudinal component of axon trajectories and reproduce the synaptic contact probabilities found in experiments (Li et al., 2007). We then used this axon growth model to generate model tadpole networks which we analyzed anatomically but not physiologically (Borisyuk et al., 2008, 2011). The present modeling is based on more detailed anatomical analysis of a larger set of filled neurons, new insights into the channel properties of spinal neurons (Winlove and Roberts, 2012), and the properties and roles of excitatory reticulospinal dIN neurons (Moult et al., 2013). It also uses a new, biologically realistic, gradient field-dependent model of axon growth.
Recordings from tadpole neurons define their properties and connections
Motor nerve recordings in immobilized tadpoles (Fig. 1D–F) show that a skin stimulus can initiate motor nerve activity (43.6 ± 11.5 ms after the stimulus: 155 trials in 31 tadpoles) matching real swimming behavior in frequency (10–25 Hz) and progressing from head to tail (Kahn et al., 1982).
Whole-cell recordings have been used to trace the chain of events from a skin stimulus to swimming. Overall, the accumulated recordings from >1500 pairs of dye-filled neurons define the properties and synaptic connections of seven neuron types forming a network that allows the tadpole to swim when touched (Fig. 1C; Soffe et al., 2009; Roberts et al., 2010; Li, 2011; Buhl et al., 2012). Crucial components of the swim network are the excitatory pacemaker dINs, which drive swimming activity. We knew that these were excited when a few sensory RB neurons were excited to spike 4.5–8 ms after local trunk skin stimulation (Soffe, 1991), but the detailed timeline for activation of the dINs by trunk stimulation was unclear. We therefore made single or paired whole-cell recordings to define the delays in sensory pathway and dIN neuron responses following trunk skin stimulation above the threshold for evoking swimming monitored by motor nerve recordings (74 tadpoles, with 17 neuron pairs). Sensory pathway neurons on the stimulated side (dlc and dla; Fig. 1C) fire shortly after the sensory RB neurons (7.6 ± 1.8 ms after skin stimulation; 16 neurons). The dINs on each side of the body then receive excitation (13.6 ± 8.1 ms, 58 neurons). This excitation leads to dIN firing on one side or the other after 36.5 ± 23.1 ms (48 neurons), and when this occurs, rhythmic swimming starts (Fig. 1E,F). In this way we have defined the input and output of the whole swim network as well as the precise patterns of neuron activity during swimming (Roberts et al., 2010).
Anatomy of seven brainstem and spinal cord neuron types in the swim network
Detailed anatomical analysis is required to allow us to model the growth of each neuron type, the crux of our developmental strategy, and the formation of the swim network in the young tadpole. The spinal cord has a simple tube-like architecture with a ventral “floor plate” formed of cuboidal cells and a dorsal “roof plate” formed by sensory RB neuron somata (Fig. 2A,B). As in all vertebrates, sensory neurons and functions are located dorsally, while motoneurons and motor functions are more ventral. At least seven types of spinal and hindbrain neuron are active during swimming (Li et al., 2001; Roberts et al., 2010). Using 276 single neuron dye fills from previous and 8 from current whole-cell recordings, details of their soma positions, dendrite dorsal and ventral extents, complete axonal trajectories, and branch points were defined by measuring a sample of each neuron type (Fig. 2C–H; Table 3). Axons were considered to be filled to their ends if they had a well labeled end-bulb. The characteristics of each neuron type were analyzed to define common patterns and generate a library of quantitative anatomical measurements. Most neuron somata lie medial to the marginal zone (MZ) where longitudinal axons make en-passant synapses with mainly radial dendrites. Most axons grow ventrally into the MZ before turning to grow longitudinally on the same side or growing though the floor plate to reach the opposite side before turning longitudinally (Fig. 2B).
The present data on the DV distribution of axons and dendrites extends our previous measurements (Li et al., 2007) in both number and details of measures. These distributions determine where synapses can form and the broad differences between different neuron types are illustrated here by their median DV positions (in micrometers relative to the ventral edge of the MZ; Figs. 2H, 3B,C). Sensory neuron axons lie in the dorsal tract, which is separated from the MZ by a barrier formed by a superficial row of dorsolateral neuron somata belonging to sensory pathway dla and dlc neurons (cf. adult dorsal column and dorsal horn; dl in Fig. 2A,B; Clarke and Roberts, 1984). Sensory pathway neurons have dendrites distributed in the dorsal MZ and the dorsal tract (median DV positions: dla, 109.9; dlc, 106.7) where they can contact sensory axons (Fig. 3C; Li et al., 2003, 2004a). The axons of all but the sensory neurons lie in the MZ (Fig. 2A,B), and their DV distributions peak in a mid-DV position (median DV positions: aIN, 45.9; dla, 45.4; dIN, 35.3; dlc, 32.6); in commissural interneurons (cINs) and motoneurons the peak is more ventral [median DV positions: cIN, 23.5; mn (motoneuron), 13.3; Fig. 3B]. The dendrites of the nonsensory pathway neurons have distributions extending through the MZ rather like the axons, but with aIN and motoneuron dendrites having a more ventral bias (median DV positions: dIN, 44.8; cIN, 42.7; mn, 34.8; aIN 31.0; Fig. 3C).
Modeling axon growth for each type of neuron
New anatomical data and their analysis provided the essential underpinning of an axon growth model. The first requirement was a two-dimensional growth environment. This consisted of a rectangular field. It was based on measures from sections of the hindbrain and rostral spinal cord (Fig. 2A) and was made by opening the CNS like a book along its dorsal midline for 2000 μm from the midbrain (Figs. 2B,H, and 3; Li et al., 2007; Borisyuk et al., 2011). Regions >2000 μm caudal to the midbrain were not considered because the narrowing of the spinal cord becomes significant and we have insufficient data on caudal neuron anatomy and physiology. Barriers to axon growth are formed by the edges of the whole CNS, the floor plate ventrally, a dorsolateral row of neuronal somata separating the marginal zone from the sensory dorsal tract (dl in Fig. 2A,B), and the roof plate of sensory RB somata dorsally. The floor plate is the only barrier routinely crossed by commissural axons. In addition to the barriers, the model growth environment was given three axon growth cues: a rostral-caudal polarity cue, and dorsal and ventral cues representing diffusive chemical gradients (Sanes et al., 2006), each with the same exponential form and slope (Fig. 3A). These DV gradients originate at the dorsal edge of the cord and within the floor plate, 5 μm from the ventral midline. The barriers, rostrocaudal (RC) polarity, DV gradients, and overall dimensions define the growth environment.
The axon growth model uses surprisingly few and simple rules to match the trajectories of the real axons of each neuron type. In essence, they are as follows. (1) Axons were assigned an initial axon position (same as soma position), an outgrowth angle, a branch distance and a branch angle, and a final length, all based on generalization from measurements (Table 3; see Materials and Methods). (2) Axons then extend in 1 μm steps, turn after each step in response to the three guidance cues (the angle depending on the modeled sensitivities of the growth cone to each gradient) with a small random element, and some branch to grow a second axon once the first has grown. An optimization process was used to find parameters for guidance cue sensitivities and the degree of random turning for axons of each neuron type giving the best match to real axon trajectories. We found that axons of each neuron type also needed a sequence of different parameter sets at successive stages of their growth. All axons have an initial orientation to longitudinal growth stage, and a main growth stage; commissural axons also have a change in their response to gradients after growing through the floor plate (Moon and Gomez, 2005). After the parameter optimization for each stage, the model generated sets of realistic axonal trajectories and distributions for each neuron type, either to the same or the opposite side of the CNS (Fig. 3B,D).
Modeling the formation of synaptic connections and the swim network
The axon growth model was then used to produce a map of connections from axons onto dendrites in the whole swim network. Defined populations of each neuron type (∼30–200) were distributed along each side of the 2D growth environment based on available evidence, for example from transcription factor immunocytochemistry (Li et al., 2001, 2004b). The axon growth model then generated the axons of the ∼1400 neurons in our 2000-μm-long model based on the simple gradients and barriers. Since understanding of spinal neuron dendrite growth is lacking (Cline and Haas, 2008; Katsuki et al., 2011), we assigned a single straight “dendrite” to each neuron at its soma position with dorsal and ventral limits generalized from measurements (Figs. 2B, 3A,C, 4B; Table 3; Li et al., 2007). Where an axon crossed a dendrite, a synapse was made with a probability of 0.46 in the marginal zone and 0.63 for sensory connections in the dorsal tract based on evidence from paired whole-cell recordings (Li et al., 2007). Each run of this process generated a unique network: a map and list of all the synapses (range: 81,822–91,045; mean number 86,655 ± 1412 SD) between all the neurons in the model network (Fig. 4A; Table 4).
Each network generated by the growth model (we will use this term as a shorthand for: axon growth and connection probability model) reflects the particular distribution of neurons, axons, and dendrites (Fig. 4). Thus, as in life, sensory RB axons in the dorsal tracts synapse with sensory pathway dlc and dla neurons, the only neurons with dendrites lying in this tract (Fig. 4C). These sensory pathway neurons in turn synapse onto the swim neurons like dINs. Since their long ascending axons accumulate rostrally, they make more synapses here and fewer caudally (Fig. 4C, dla and dlc; Wolf et al., 2009). The excitatory dINs (∼120 per side), which drive swimming, make synapses with each other, with inhibitory cIN and aIN neurons, and with motoneurons. The reciprocal inhibitory cIN neurons make synapses with all except sensory RB neurons on the opposite side (for summary see Fig. 1C; Table 4)
Can modeled networks produce realistic, whole animal behavior?
Crucially, we can validate and test the adequacy of each “grown” network because we know the real network response to skin stimulation is sustained, coordinated swimming. The generated swim networks were therefore mapped onto a physiological model to ask if a brief stimulus to some of the sensory neurons innervating the skin triggers sustained swimming. Neurons in the physiological model were single-compartment and, for simplicity, most of them had the same, typical, Hodgkin–Huxley-type membrane channels as motoneurons in our previous study (Dale, 1995; Sautois et al., 2007). The model dIN neurons alone were given special properties based on recent experimental evidence from whole-cell paired and voltage-clamp recordings and modeling using a population of electrically coupled multicompartment dINs (Roberts et al., 2010; Li, 2011; Winlove and Roberts, 2012). The majority of voltage-gated current is carried through one type of sodium channel (Na), fast and slow potassium channels (Kf, Ks), and a high-voltage-activated calcium channel (Ca). After tests in the multicompartment model, the channel properties and weak electrical coupling were mapped onto single-compartment dINs in a network model. Injecting current into a single dIN showed that the new model neurons responded with single spikes but could show postinhibitory rebound firing as seen in experiments (Fig. 5A,B; Li et al., 2006). However, we could also show that the electrically coupled population of dINs would show rhythmic pacemaker-like firing in the swimming frequency range during the slow activation of an NMDA conductance (Fig. 5C; Li et al., 2010). In the swimming networks we build, rhythmic firing is based on feedback glutamate excitation activating NMDA and AMPA receptors. It depends on voltage-gated Na+, Ca2+, and fast and slow K+ channels and postinhibitory rebound firing. We used three types of conductance-based synapses: glutamatergic AMPA and NMDA, and glycinergic. Synaptic time course and strengths were based on evidence from paired recordings (Dale, 1995; Sautois et al., 2007) with only minor adjustments to connections strengths which are defined in Materials and Methods.
When a network was mapped onto a physiological model incorporating the special dIN properties, a single spike in two sensory RB neurons (modeling touch to the skin) initiated the alternating motoneuron firing we define as “swimming” (Fig. 6A–C). Stimulation made sensory neurons (yellow/black) fire and recruited sensory pathway neurons (pink and red). These amplified the sensory excitation and distributed it to both sides. Swimming started first reliably on the nonstimulated side, provided that sensory pathway neuron synapses were asymmetric and dlc neurons made stronger synapses on that side (Zhao et al., 1998). Alternating swimming activity occurred reliably in most networks tested (99 of 100 tests), often preceded by several cycles of synchronous left-right firing at double the swimming frequency (82 of 99 tests; range 1–7 cycles, the mode is 1 cycle). The swimming frequency (17.8 ± 0.55 Hz) was within the normal range (10–25 Hz) and left/right motoneuron activity was strictly alternating (phase 0.5 ± 0.01, n = 99) as in the immobilized tadpole (Kahn et al., 1982). If sensory pathway dla neurons projecting to the stimulated side made stronger synapses, motoneuron activity started on the stimulated side (Fig. 6D; 12 of 12), but there was little effect on average swimming frequency (18.2 ± 0.4 Hz). Without such asymmetry the network still swam reliably (11 of 12 tests) with an average frequency 18.1 ± 0.5 Hz, but swimming was preceded by more synchronous motoneuron activity on both sides of the body after stimulation (11 of 12 tests; range 2–7 cycles, the mode is 7 cycles). Head-to-tail progression of motoneuron activity was not clear in the model. However, it was also not found in the rostral 2 mm in vivo and only becomes clear over longer regions of the spinal cord (Tunstall and Roberts, 1991). The importance of the pacemaker properties of dINs was emphasized by the observation that if dINs were given the same properties as other neurons, a reflex occurred, but rhythmic swimming failed (0 of 12; Fig. 6E). In vivo, the initiation pathway is probably more complex because the delay to the first motoneuron spike is longer (43.6 ± 11.5 ms vs 19.6 ± 0.5 ms in the model) and can be on either side (50:50 in 155 trials in 31 tadpoles). Synchrony following stimulation is rare in vivo. However, these results show that complex model networks built by generalizing from small biological datasets and assembled following remarkably simple rules can produce reliable coordinated motor activity like swimming. The patterns of activity shown by each type of neuron are also similar to those seen in whole-cell recordings (Figs. 1E,F, 6B).
Experiments to find the essential features of model network assembly
The physiological model was then used to test the effect of easing anatomical constraints on network self-assembly. Swimming remained reliable but became faster (23.2 Hz ± 0.6, n = 12) when synapse formation probability was increased to 1 from 0.46 and 0.63. Swimming was still reliable when synapse probability was reduced by 12% (12 of 12). After a 25% reduction, most networks swam (9 of 12), whereas the rest showed rhythmic activity on one side of the body only. Scaling the probability down by 25% again (yielding a 44% reduction from the original value) eliminated swimming in all networks (0 of 12). We then tested whether the details of the axon or dendrite distributions of the neurons active during swimming (dIN, cIN, aIN, mn) in the marginal zone were significant for swimming. Swimming was still reliable when these neurons were all given the same axon DV distributions as sensory pathway dlas (at 16.2 ± 0.5 Hz; 12 of 12), or the dendrite DV distributions of inhibitory cINs (at 17.9 ± 1.0 Hz; 12 of 12; Fig. 7A). The next test was to make all CPG neuron dendrites span the whole dorsoventral extent of the marginal zone so they could be contacted by all axons except those of sensory neurons. This test effectively discards all the details of the axon trajectories and dendritic arborizations of each neuron type. Remarkably, these networks produced reliable swimming provided that synapse formation probability for CPG neurons was reduced (to 0.25) so that total synapse number was conserved (swimming at 13.9 ± 0.3 Hz; 12 of 12; Fig. 7B,C).
If the detailed DV positions of axons and dendrites are not critical to swimming network function, we can extend our model to more caudal regions of the spinal cord where anatomical data are scant. We therefore extended the model from 2000 to 3500 μm so it had 1980 neurons. Stimulation produced swimming as previously shown (at 17.5 ± 0.4 Hz; 12 of 12; Fig. 7D,E) but a head-to-tail delay of 4.4 ms mm−1 was now clear between the earliest motoneuron spikes (at ∼1000 μm) and the latest (at ∼2500 μm). This is within the range of 2–5 ms mm−1 found experimentally (Tunstall and Roberts, 1991).
The final test was to make the dendrites of all CPG neurons in the 2000 μm network reach into the dorsal tract so they could be contacted and excited by sensory RB axons. Sensory stimulation now led to short-latency firing in all types of neurons on the stimulated side, which was never seen in recordings. Swimming always followed (at 19.1 ± 0.5 Hz; 12 of 12) but was preceded by more synchrony (mean 3.8 cycles) than in the control case.
Discussion
We show that it is possible to use a model of neuronal development to generate the large-scale anatomical pattern of neurons, axons, and synaptic connections forming the core of the vertebrate nervous system. This also leads to the conclusion that the first functional networks in the vertebrate brainstem and spinal cord may develop using surprisingly simple rules. This suggests that complex and large networks can assemble where connections are made without recognition of “correct” target neurons. Axons may only need to distinguish dendrites from glia and axons. They can then synapse with any dendrite they contact (with a certain probability) so long as they grow into an appropriate region and in a broad direction along the nervous system (toward the head or tail; on the same or the opposite side). This important finding has implications for the longstanding debate on how neuronal connections are made (Sperry, 1963; Zipursky and Sanes, 2010). First, the circuits in developing larval vertebrates are not simply the precursors of more complex networks that only become effective once their connectivity has been refined by additional developmental processes like detailed recognition. They have to function immediately in their own right. As soon as they hatch, fish and frog larvae need to swim to avoid predation, so the networks controlling their first behavior have to work properly. Second, if the formation of spinal circuits relies on simple processes, then there is hope for regeneration of spinal circuitry after injury if injected precursor cells only require some simple guidance to navigate the injured environment and integrate into existing circuitry (Bonner et al., 2011; Tuszynski and Steward, 2012).
Simple structural features may be critical for network formation. In the adult spinal cord, sensory axons enter the dorsal horn and dorsal column and are separated from more ventral neurons organizing action. This fundamental segregation of sensory and motor functions is also a feature of the developing frog spinal cord but has received little previous attention. The superficial somata of sensory pathway neurons (in a position corresponding to the adult dorsal horn) separate skin sensory axons in the dorsal tract (Nordlander, 1984, 1987; in a position corresponding to the spinal dorsal column) from the more ventral axons and dendrites of the other neurons. The consequence of this barrier is that sensory axons only make synapses with sensory pathway neurons, because these are the only neurons with dendrites in the dorsal tract. While this segregation of sensory functions is important, our manipulations of model dendrite lengths show that detailed specification of the dorsoventral distribution of axons and dendrite lengths in our models may actually not be necessary for swimming (Fig. 7B,C). However, the broad pattern of longitudinal neuron distributions, axon projections (ipsilateral or contralateral, ascending and/or descending), and axon lengths are important because they determine the distribution of synapses along the body (Wolf et al., 2009). Our modeling points to these basic features and transmitter phenotypes controlled by transcription factor expression (Lewis, 2006) being necessary in the neuron specification process which can lead to early network formation.
Our study establishes that simple axon guidance mechanisms may be sufficient for the assembly of networks producing the first behavior of developing vertebrates. In fact, the real mechanisms directing axon growth may not be the ones we have used to generate our functional networks (Forbes et al., 2012). In a number of systems there is detailed information on the mechanisms guiding growing axons to their target areas (Dickson, 2002; Chilton, 2006; Zou and Lyuksyutova, 2007; Kastanenka and Landmesser, 2010) and the role of morphogen gradients in this process (Arber, 2012). To allow highly specific synaptic connections, the initial broad mapping of axons is often refined by other mechanisms, sometimes dependent on the timing of neuron spiking activity (Kastanenka and Landmesser, 2010). These processes lead to remarkable examples of precision axon-mapping like retinal ganglion cell projections to the brain (McLaughlin and O'Leary, 2005), motoneurons to specific limb muscles, and muscle stretch receptor afferents to motoneurons (Pecho-Vrieseling et al., 2009). For simplicity, we chose chemical gradients as classic developmental cues (Sperry, 1963) and allowed each neuron to grow its axons independently. This is not what happens in life, where neurons differentiate in a defined dorsoventral and head-to-tail sequence (Lewis, 2006; Chédotal and Richards, 2010). As a result, the axons of more precocious neurons can be followed by fasciculation of those that develop later. The lack of precision in axon trajectories implied by our modeling still leaves many details to be explained, as is apparent from the list of parameters we needed to assign: initial outgrowth angle, orientation to longitudinal growth, distance to axon branching and its angle, axon length (Table 1). With the exception of axon length, control of all these variables occurs physically quite close to the soma. Does this suggest that later control of axon growth is minimal? Is dendrite growth controlled precisely or could it be a simple response to the presence of axons ready to make synapses (Cline and Haas, 2008; Arikkath, 2012)? Details of the real processes now need to be investigated experimentally in the light of our results.
Is there a role for subtlety and complexity when the first functioning networks develop in the vertebrate brain? Even in early larval swimming networks there is evidence that homeostatic mechanisms regulate synaptic strengths (Borodinsky et al., 2004), switches occur in the mechanisms controlling axon growth (Moon and Gomez, 2005), the order of neuron differentiation may affect recruitment during swimming (Koyama et al., 2011), and detailed physiological properties of neurons are critical for swimming (Li, 2011). It is highly likely that refinement of connections occurs later in development, as in most systems studied (Sanes et al., 2006), but our results suggest that it is not required for early network function. Specific target recognition may not be necessary until limbs and eyes develop. If simple rules lay out the first scaffold of axons and connections in the vertebrate hindbrain and spinal cord, then the first steps in the construction of networks in other brain regions, even the cerebral cortex (Hill et al., 2012), could be based on similar principles.
What is the significance of this study for understanding of neuronal rhythm generation networks? When the tadpole CNS network was mapped onto a physiological network model where the main excitatory neurons have pacemaker properties, swimming motor output was a remarkably stable response to sensory stimulation. This is what a very young animal may need to survive (at this stage of development, the hatchling tadpole does not appear to have a fast C-start response). The stability almost certainly results from the incorporation of recently described switchable pacemaker properties into models of the dINs which drive swimming (Li et al., 2010). The dINs in the hindbrain may be homologues of reticulospinal neurons in adults. Neurons with related pacemaker properties have been found in many central pattern generators underlying rhythmic activity, where they act in concert with network-based rhythm generation mechanisms (Marder et al., 2005; Huss et al., 2008; Feldman et al., 2013; Moult et al., 2013). In the tadpole, many issues still remain. What is the cellular basis of dIN pacemaker properties? Why does network activity start at shorter latencies in the model networks than in the tadpole? Is there any biological sense in the periods of synchronous activity on both sides of the body (Fig. 6C,D) which are also seen in physiological recordings and other models of reciprocally inhibitory networks (Wang and Rinzel, 1992)?
In conclusion, this study has used novel methods to recreate a significant part of the connectivity of the brainstem and spinal cord of a whole vertebrate by generalization from relatively small anatomical datasets defining the features of the main constituent neuron types. Such data, critically including axon trajectories, are far from ideal or complete but are just not available for other animals beyond Caenorhabditis elegans (White et al., 1986), and even here, gaps need to be filled by extrapolation (Haspel and O'Donovan, 2011). Using a simple growth model, we have been able to exploit limited data to generate a biologically realistic baseline network for a significant part of the newly hatched Xenopus tadpole CNS. The details of the distribution of the ∼90,000 synapses between ∼1400 neurons (Table 3; Fig. 4) are now available for experimental and anatomical testing, and we have already begun to ask what is necessary to ensure proper network function (Fig. 7). In the tadpole we have been able to combine generated networks with detailed physiological data to define the structure and function of a network driving behavior. This process gives insights into the formation as well as the tolerance to variability and operation principles of such networks. It is only possible because we have detailed information on the neuron types, their anatomy, properties, and synaptic connections. Such details are gradually emerging for the related networks in other vertebrates (Fetcho and McLean, 2010) including mammals (Goulding and Pfaff, 2005; Grillner and Jessell, 2009).
Footnotes
This work was funded by the UK Biology and Biotechnology Research Council (BB/G006652/1, BB/G006369/1). Studentships to support M.H. and R.M.-H. came from EPSRC/MRC/BBSRC-funded Doctoral Training Centre for Neuroinformatics at Edinburgh and School of Computing and Mathematics, University of Plymouth. R.B. is on leave from Institute of Mathematical Problems in Biology, Russian Academy of Sciences, Pushchino, Moscow Region, 142290, Russia. We thank Antony Dodd, James Hodge, Mark Holderied, Conor Houghton, Wen-Chang Li, Robert Meech, Nick Roberts, Shelby Temple, and David Wilshaw for helpful comments on the manuscript, and Debbie Carter, Tim Colborn, and Lia Vincent-Gilmour for technical assistance.
The authors declare no competing financial interests.
- Correspondence should be addressed to Alan Roberts, School of Biological Sciences, University of Bristol, Woodland Road, Bristol, BS8 1UG, UK. a.roberts{at}bristol.ac.uk