Structural and functional properties of a probabilistic model of neuronal connectivity in a simple locomotor network
 Cited 5
 Views 1,274
 Annotations
Abstract
Although, in most animals, brain connectivity varies between individuals, behaviour is often similar across a species. What fundamental structural properties are shared across individual networks that define this behaviour? We describe a probabilistic model of connectivity in the hatchling Xenopus tadpole spinal cord which, when combined with a spiking model, reliably produces rhythmic activity corresponding to swimming. The probabilistic model allows calculation of structural characteristics that reflect common network properties, independent of individual network realisations. We use the structural characteristics to study examples of neuronal dynamics, in the complete network and various subnetworks, and this allows us to explain the basis for key experimental findings, and make predictions for experiments. We also study how structural and functional features differ between detailed anatomical connectomes and those generated by our new, simpler, model (metamodel).
https://doi.org/10.7554/eLife.33281.001Introduction
Information processing in the brain is based on communication between spiking neurons that are embedded in a network of synaptic connections. Clarifying the interplay between network connectivity and functionality is a key part of understanding how the brain generates functional behaviours (Sporns et al., 2005; Marder and Calabrese, 1996). Studying this relationship is difficult because nervous system connectivity usually varies considerably between individuals. Despite this variation each individual behaves in approximately the same way, especially in the case of simple animals. This commonality of behaviour suggests that there are some fundamental organisational principles that underlie the structure of a species’ nervous system. How can we identify these fundamental properties that are shared across individuals and allow the nervous system to function correctly?
In this paper, we attempt to answer this question in the case of the hatchling Xenopus tadpole. Whole cell recordings and anatomical measurements of neurons, combined with computational modelling, have uncovered many important details regarding the neuronal network that controls swimming in hatchling tadpoles (Roberts et al., 2010; Roberts et al., 2014). We have previously shown how modelling of the neuronal connectivity in the tadpole spinal cord and caudal hindbrain is possible through a ‘developmental’ approach, whereby connections between neurons are not prescribed but appear as a result of the intersection between (simulated) growing axons with dendrites (Borisyuk et al., 2014). This anatomical model mimics the realistic growth of axons in the spinal cord. Following biological realism, the axon growth is guided by the concentration of chemical gradients in the spinal cord. The properties of such gradients are controlled by model parameters that have been optimized to produce the same statistical characteristics as real measurements. Other model specifications (including soma positions and dendritic extents) are assigned from the distributions of experimental data and from general biological knowledge. The model includes several stochastic components (Borisyuk et al., 2014; Roberts et al., 2014); therefore, each model simulation generates a different pattern of connectivity (‘connectome’). The connectivity can be mapped onto a functional model composed of spiking units of HodgkinHuxley type, with parameters chosen to match known tadpole electrophysiology (Sautois et al., 2007). The resulting functional model reliably produces activity patterns like those seen during real swimming (Roberts et al., 2014). It is important to note that the anatomical model provides a way of generating many different connectomes, such that the random variation observed between generated connectomes has the same statistical properties as measurements taken from different individual animals. Here, we set out to reveal the fundamental features of the neuronal connectivity that underlie the ability of the swim network to function robustly.
We describe a new probabilistic model of connectivity, which is generalised from a large number of connectomes generated by the anatomical model. This probabilistic model is a matrix that specifies the probability of connection between each pair of neurons. Being derived from multiple biologically realistic (anatomical) connectomes, the probabilistic model reflects the anatomical structure of the biological system. An important advantage of the probabilistic model is that it is simple enough that we can analyse the properties of the model itself, rather than individual connectome realisations. We use the probabilistic model to calculate structural properties of the tadpole network. These results are general, and therefore should reflect the fundamental organisational principles that we aim to uncover here.
Graph and network theories (Rubinov and Sporns, 2010) are increasingly used to study connectivity of different neuronal networks: C. elegans (Varshney et al., 2011; Kaiser and Hilgetag, 2006), zebrafish (Stobb et al., 2012), cat, rat and macaque cortical structures (Sporns et al., 2007; Sporns and Zwi, 2004; Humphries et al., 2006). For example, it was shown that the C. elegans connectome is heterogeneous and has a hub structure (Towlson et al., 2013): most neurons have a low number of connections but there are several highly connected ‘hubs’. Hubs have been identified in many brain networks and they are likely to be formed at an early stage of development (Varier and Kaiser, 2011). However, not all brain circuits have hubs; for example, they have not been found in the rat reticular formation (Humphries et al., 2006) or zebrafish nervous system (Stobb et al., 2012).
Using the probabilistic model, we estimate the heterogeneity (Hu and Wang, 2008) and connection degree distributions (Barabasi and Albert, 1999; Sporns et al., 2007; Varshney et al., 2011) of the tadpole’s spinal cord network. We found that the generalised tadpole network is not scalefree and that hubs do not exist; therefore in this respect the generalised tadpole network differs from the C. elegans connectome.
A second potential advantage of the probabilistic model is that it can be used to easily generate connectome realisations by sampling from the probability matrix, without detailed simulation of neuronal growth. This enhances its potential value as a tool for studying the functional properties of the network when combined with an appropriate physiological model. Multiple functional simulations of probabilistic connectomes demonstrated a reliable pattern of rhythmic activity, qualitatively like tadpole swimming and as seen in previous modelling (Roberts et al., 2014). Thus, the generalised probabilistic model shares structural and functional properties with the real biological object. However, quantitative differences showed that caution is required to avoid pitfalls when employing the probabilistic approach to study real biological activity.
Specifically, we found that the variance of the number of incoming connections (indegree) or outgoing connections (outdegree) of each neuron is higher in anatomical rather than a probabilistic connectomes. As a result of this finding, we observed that the period of the rhythm was longer in probabilistic connectomes. We can explain why the generalisation process affects the swimming period, and show how it is possible to accurately predict the period of swimming using only structural properties of the connectome. We then show how, by making suitable parameter adjustments, we can match the functionality of the probabilistic connectomes to that of the animal and anatomical connectomes. This makes it possible to use the probabilistic approach as a tool for studying real biological activity as well as structural properties of networks.
Despite the differences between anatomical and probabilistic models, we demonstrate several important advantages of using the probabilistic model in comparison to the anatomical one. For example, we could predict the position of commissural interneurons (cINs) that are active during swimming, which can be difficult to explain by the anatomical model. Specifically, our simulations show that cINs in rostral positions are less likely to fire reliably than those in caudal positions. Moreover, the probabilistic model allowed us to easily design new computational experiments that helped to clarify the following experimental findings.
By studying the connectivity of CPG neurons specifically, we show that the minimal swimming subnetwork includes neurons of two types: inhibitory commissural interneurons (cINs) and excitatory descending interneurons (dINs). Similar to experiments with the surgically isolated half semiCNS (Soffe, 1989), we found that the network of interconnected dINs on one body side could still generate rhythmic activity even without commissarial inhibition. It is known from experimental measurements that some dINs have both descending and ascending axons (Roberts et al., 2010). Our simulations of the model without ascending dIN axons show that the ascending connections play a key role in swimming and their deletion leads to pathological activity.
To summarize, in this paper we design a simple probabilistic model (metamodel) which reflects some structural features of anatomical connectomes. We also show that it can be used to study how these features relate to real behaviour by making suitable adjustments in synaptic strengths. We consider this investigation of the tadpole spinal cord as an important example of a technique that can be widely applied to study the nervous system of other animals.
Materials and methods
Derivation of the probabilistic connectivity model
Request a detailed protocolThe probabilistic connectivity model is derived from multiple connectomes generated by our existing anatomical model: a developmentally inspired model which is biologically realistic and incorporates a large number of biological measurements (Borisyuk et al., 2014; Li et al., 2007a; Borisyuk et al., 2011). The anatomical model simulates axon growth guided by chemical gradients, with model parameters that are chosen by fitting the generated axons to experimental measurements. As the growing axons intersect dendrites, which are allocated along the body according to experimental measurements, synapses form and make connections between neurons.
Here, we explain some details of the anatomical model that are important for understanding the new probabilistic model (for more details about the anatomical model, see Appendix 1). The anatomical model includes N = 1382 neurons of the seven types known to generate the swimming response. The network is divided between neurons in the sensory pathway (RB, dlc and dla), CPG neurons (dIN, cIN and aIN), and output motor neurons (mn). Sensory pathway neurons deliver sensory stimulation to CPG neurons. CPG neurons are responsible for the generation and maintenance of the swimming activity pattern. Motor neurons (mn) deliver CPG output to muscles and generate locomotion (Figure 1). The model is simplified by fixing the number of cells for each neural type, with neurons of each type equally divided between the left and right body sides. Simulation of the anatomical model results in a network with approximately 83,000 synapses on average. For a full description of the anatomical model and its implementation, see (Borisyuk et al., 2014; Roberts et al., 2014).
Importantly, the anatomical model includes stochastic components, so repeatedly running the model produces different connectomes with different numbers of connections and connection distributions. In particular, rostrocaudal coordinates of neurons can vary between connectomes. However, since the number of neurons of each type is kept constant it is possible to find a onetoone correspondence between any two generated connectomes. First, we ordered the cell types (RB, dla, dlc, aIN, cIN, dIN, mn) and second, for each cell type we ordered neurons of that type according their longitudinal position (or the rostrocaudal (RC) coordinate) in ascending order from head to tail. For example, in any connectome neuron #1 is the most rostral RB neuron on the left side of the body, while neuron #62 is the most caudal leftside RB; neurons #63–126 are the rightside RB neurons; neurons #127–146 and neurons #147–174 are the dla neurons on the left and right sides respectively, etc. This ordering of cells is universal and does not depend on a particular connectome; therefore, we can enumerate all neurons in a universal way, providing a onetoone correspondence between generated connectomes.
To define the probabilistic model we used the universal enumeration of neurons and considered the matrix of probabilities $P$ where ${p}_{i,j}$ is the probability that there is a synaptic connection from neuron $i$ to neuron $j,i=1,2,\dots ,N;j=\mathrm{1,2},\dots ,N.$ Here N = 1382 is the total number of neurons. We defined the random Bernoulli variable ${X}_{ij}\in \left\{\mathrm{0,1}\right\};$ where ${X}_{ij}=1$means that there is a directed connection from $i$ to $j$ and the probability $\mathrm{P}\mathrm{r}\{{X}_{ij}=1\}={p}_{i,j}$. To calculate an estimate of this probability ($\hat{p}}_{ij$), we generated $K=1000$ connectomes and calculated the frequency of appearance of this directed connection: $\hat{p}}_{ij}=\frac{M}{K$, where $M$ is the number of connectomes with a connection from neuron $i$ to neuron $j.$ The RCcoordinate of each neuron is defined by the averaging the RCcoordinates across the $K$ generated connectomes.
The central limit theorem provides the error estimation of each entry of the probability matrix ${p}_{ij}$: the length of the binomial confidence interval with 95% confidence is given by $\text{}{e}_{ij}\approx 2\frac{1.96}{\sqrt{M}}\sqrt{{\hat{p}}_{ij}\left(1{\hat{p}}_{ij}\right)}$. The maximum of this error’s estimate corresponds to ${\hat{p}}_{ij}=0.5$, therefore, $\underset{i,j}{\mathrm{max}}\left({e}_{ij}\right)\approx 0.06.$
The probabilities of directed connections between all neurons of the swimming network are shown in Figure 2A. All probabilities are between ${\hat{p}}_{ij}=0$ (no connections) and ${\hat{p}}_{ij}=0.69$. To visualize these probabilities we use a greyscale, where black color corresponds to ${\hat{p}}_{ij}=0$ and bright pixels to high probabilities. Note: here and below we use the same notation ${p}_{i,j}$ for the probabilities and their estimates.
As an example, Figure 2B shows the submatrix corresponding to aINaIN connections. There is a black diagonal line which results from the fact that neurons cannot make connections with themselves. In fact, similar almostdiagonal lines can be seen in all of the other submatrices due to a feature of the growth model that prevents neurons contacting very nearby neurons. Close to the diagonal line in Figure 2B the shading is very bright, but this fades to black away from the diagonal. This results from the fact that the probability of two neurons being connected decreases with the distance between them. For aINs, the shading is brighter below the diagonal line, which reflects the fact that their axons are mainly in the ascending direction, making a given aIN more likely to contact aINs that are located more rostrally. While the aINaIN example is relatively simple to understand, neurons with more complicated growth patterns have submatrices with more complex structure – for example in the case of dINdIN connections.
The matrix $P$ can be used to generate a specific adjacency matrix of directed connections (connectome) $A=\left({a}_{ij}\right)$ where ${a}_{ij}\in \left\{\mathrm{0,1}\right\}$ and ${a}_{ij}=1$ indicates existence of the connection from neuron $i$ to neuron $j$. This matrix A is a particular realization of independent Bernoulli variables. We then used these specific adjacency matrices (‘probabilistic connectomes’) to explore their functional properties by mapping the connectomes onto our functional model to study the spiking activity in the swim network in response to stimulation.
Functional model of spiking activity
Request a detailed protocolTo investigate the relationship between the network’s structure and functionality it was necessary to simulate the spiking activity using the connectomes generated by the probabilistic model. We produced specific adjacency matrices from the probability matrix and used them in a functional model to simulate responses to stimulation and study spiking activity patterns. The functional model included conductancebased singlecompartment neurons of HodgkinHuxley type with synaptic and axonal delays. In addition to the chemical synapses that are generated by the anatomical or probabilistic models, we also included the effects of electrical coupling (gap junctions) between dINs that are in close proximity to each other. We follow previous experimental (Li et al., 2009) and modelling (Hull et al., 2015) studies that have suggested that these electrical connections are an important functional property of the dIN network. Synaptic strengths, membrane channel conductances and neuron capacitances were all based on experimental results and then randomised according to a Gaussian distribution. A complete description of the functional model is given in (Roberts et al., 2014).
Simulations were performed using NEURON 7.3 (Carnivale and Hines, 2006) (RRID:SCR_005393) with a fixed timestep of 0.01 ms.
Details about the functional model and parameter values are given in the Appendix 2. The code for the anatomical, probabilistic and functional models, and the code for reproducing all the figures in this manuscript are available in Model DB from https://senselab.med.yale.edu/ModelDB/enterCode.cshtml?model=238332.
Results
In and outdegrees derived from the probabilistic model
One way the structure of a network can be measured is by calculating the number of incoming and outgoing connections each element in the network has. In this section, we use the probabilistic model to calculate the mathematical expectation of the incoming connection number (indegree) and outgoing connection number (outdegree) (Bullmore and Sporns, 2009) for the whole network and for different subnetworks.
Based on the assumption that the probability matrix $\text{}P=\left({\hat{p}}_{ij}\right)$ consists of independent Bernoulli random variables, the mathematical expectation of the indegree ${I}_{j}$ and the outdegree ${O}_{j}$ for neuron $j$ are given by the following formulas:
These formulas follow from the fact that the random variables ${I}_{j}$ and ${O}_{j}$ have the Poisson binomial distribution (Sprott, 1958). Similarly, the formulas for the standard deviation of these random variables are the following:
Figure 3A shows the mathematical expectation and standard deviation (calculated using formula (1) and (2)) of the indegree (upper panel) and outdegree (lower panel) for each neuron on each body side according the RCcoordinate (the equivalent figure for the right side is very similar and omitted here). For each cell type, the shape of the in and outdegree distribution is very specific and depends on the soma position. For example, motor neurons (green), have high indegree and very low outdegree. Almost all shapes are unimodal with a skewed position of the maximum. This is a consequence of the interplay between primary and secondary axons in the developmental model and the RCcoordinate distributions of their somata. The absence of descending axons for dla neurons (pink) and their ‘parallel’ pattern of growth leads to almost linear increase of their outdegree. A similar explanation applies to the linearly decreasing shape of in and outdegrees for dINs near the tail, with RCcoordinates more than 1400 $\mu m$, which have only descending axons. Interestingly, aINs have high in and outdegrees, suggesting that, on grounds of connectivity, they could play a significant role in the network activity; however, both experiments and simulations of the functional model revealed that aINs are rarely active during swimming. This emphasises the key importance of considering both structural and functional properties in network activity.
Figure 3B shows estimates of the mean and the standard deviation for the connectomes which are generated by the anatomical model; we numerate them by index $\alpha (\alpha =\mathrm{1,2},\dots ,1000)$. For each neuron $i$ of generated connectome $\alpha $, we consider samples of in and outdegrees:${I}_{i}^{\alpha},\text{}{O}_{i}^{\alpha}\left(i=1,2,\dots ,\text{}N=1382\right)$. We use these samples to calculate the estimates of the mean and the standard deviation. Obviously, the average in and out degrees for anatomical connectomes (black dots in Figure 3A) are exactly the same as the mathematical expectations of the probabilistic model shown by black dots in Figure 3B. However, the estimates of standard deviation for anatomical model are significantly larger for many neurons. In case of dla and dlc cell types the standard deviation of incoming connections is similar for both anatomical and probabilistic models. The reason is that the neurons of these cell types receive connections from sensory rb neurons and the number of incoming connections to dla and dlc neurons has a very low variability.
The independence of in and outdegrees when the whole network is considered together is characterised by a small value of the correlation coefficient $r=0.07$. However, some subnets showed strong dependencies. Scatter plots (pairs $\left({I}_{j},{O}_{j}\right),j=\mathrm{1,2},\dots ,N,$ where N is the number of pairs) for cINs and dINs (Figure 3C) showed the linear dependence for dINs ($r=0.99$). Other neuron types, for example cINs, show some more complicated dependence.
A key structural property of a network is whether or not it is scalefree. A scalefree network contains some ‘hub’ nodes with large numbers of connections in comparison to other nodes, and is particularly robust to removal of random nodes (Barabasi and Albert, 1999). In line with the standard approaches used for analysing scalefree structures, we calculated the distributions of in and outdegrees. We found that all these distributions are localized around the mean and they have no tail (for this reason, we do not show these distributions here). Therefore, all these networks are not scalefree. One way in which a network can be categorised as scalefree or not is by quantifying the heterogeneity of its nodes’ in or outdegrees. We calculated the socalled heterogeneity index $H$ (Hu and Wang, 2008) for in and outdistributions to estimate the variability of in and outdegrees. We compute this index to confirm that it is less than the threshold for scalefree networks (Hu and Wang, 2008). The heterogeneity is given by the following formula:
Here ${d}_{i}$ is either in or outdegree of neuron $i$, $\stackrel{}{d}=\sum _{i}{d}_{i}$ is the average degree (either in or out), and N is the number of neurons. Note, we calculate the heterogeneity index using the probabilistic model without considering any particular connectome.
Figure 4A and B show the value of $H$ for each celltype and for in and outdegrees respectively. A standard approach for determining whether or not a network is scalefree is to compare its heterogeneity index with that of a known scalefree random network. It is known that random scalefree networks with power $2\le \alpha \le 3$ have $H\ge 0.3$ (Hu and Wang, 2008). In contrast, all heterogeneity indices that were calculated (Figure 4) were relatively low (for indegrees: $H<0.2$, for outdegree: $H<0.3$); therefore, all degree distributions in the probabilistic model were rather homogeneous. From this we concluded that for each cell type the network is not scalefree, and therefore does not contain hub neurons. Thus, the connectivity of the whole tadpole spinal cord network appears organized in such a way that there are no hubs.
Another standard approach for detecting the heterogeneity and the presence of hubs is to analyse the decay of the tail of the degree distributions and compare the rate of decay with various standard functions: typically the Power Law, Exponential or Weibull distributions (Clauset et al., 2009). As previously said, in case of the probabilistic model the in and out degree histograms have no tail and therefore they differ from each of these distributions. Thus, this method for heterogeneity estimation is not applicable.
Functional properties of the model: reliable swimming
The next stage was to investigate the spiking activity of connectomes to see whether they behaved like those generated anatomically (Borisyuk et al., 2014) and as described behaviourally (Roberts et al., 2014). This was necessary to evaluate whether the probabilistic approach provided a useful tool for exploring biological function.
To investigate the spiking activity of connectomes generated by the probabilistic model, we mapped them onto a functional model composed of single compartment HodgkinHuxley type neurons, following the approach described in Roberts et al. (2014). To simulate the basic experiment where brief stimulation of the trunk skin initiates swimming in the tadpole, we excited two adjacent sensory RB neurons on one side of the body at a randomly selected RC position. The RB activity propagates along their own axons and then in the sensory pathway (via dla and dlc neurons) to deliver excitation to CPG neurons on both sides of the body. These CPG neurons (cIN, dIN, aIN) generate a pattern of rhythmic spiking alternating between the left and right body sides suitable to drive swimming movements. We repeated this experiment 100 times using different generated adjacency matrices. We found that in all simulations the functional model produced a swimminglike pattern where: firing was rhythmic; neurons that were active fired once per cycle; firing alternated between the two sides; and firing on each cycle was most delayed towards the tail.
However, although connectomes from both the anatomical and probabilistic model produced qualitatively similar swimming activity, the probabilistic model produced a rhythm with significantly longer cycle periods (68.6 ± 0.8 ms (mean ±SD), range from 65 to 70 ms) than the anatomical connectomes (58 ± 1.8 ms), as shown in Figure 5A. We investigated the underlying cause of this difference, and in doing so gained an insight into how the structure of the network affects swimming period, a key characteristic of the system’s behaviour.
What determines the period of one swimming cycle? A swimming cycle starts when dINs on one side of the body spike. These excite cINs on the same side, which then spike and inhibit dINs on the opposite side, leading to delayed spiking of dINs on the opposite side through postinhibitory rebound (PIR). Thus, the swimming period can be approximated as $T=2\left({\mathrm{\Delta}}_{DC}+{\mathrm{\Delta}}_{CD}\right)$, where ${\mathrm{\Delta}}_{DC}$ is the delay between spiking of dINs and the subsequent spiking of the ipsilateral cINs they excite, and ${\mathrm{\Delta}}_{CD}$ is the delay between spiking of cINs and the subsequent PIR spiking of the contralateral dINs they inhibit (Figure 5B). Both ${\mathrm{\Delta}}_{DC}$ and, particularly, ${\mathrm{\Delta}}_{CD}$ were significantly larger with the probabilistic connectome (anatomical model: ${\mathrm{\Delta}}_{DC}=5.3ms\pm 0.4,{\mathrm{\Delta}}_{CD}=23.7ms\pm 0.9,N=100$; probabilistic model: ${\mathrm{\Delta}}_{DC}=6.2ms\pm 0.3,{\mathrm{\Delta}}_{CD}=28.2ms\pm 0.4,N=100$). Together these two differences account for the overall slower swimming rhythm seen with the probabilistic model, and this is largely as a result of the increased time it takes for dINs to fire PIR spikes in response to contralateral cIN input.
What, then, determines the delay between cIN spikes and contralateral dIN rebound spiking? During swimming dINs are held depolarized by summation of NMDAreceptormediated excitation from other dINs, and in this state inhibition from cINs can result in delayed dIN spiking as a result of PIR. Intuitively, and from past investigations, we know that this spiking delay depends on the relative strength of inhibitory and excitatory input from cINs and other dINs, respectively. We characterised the relative strength of inhibition and excitation for a given connectome by calculating the average indegree from cINs and from other dINs. Any cINs that received fewer than 13 connections from dINs were excluded from this calculation, since, as we shall demonstrate, such cINs are likely to be inactive. We used a linear regression model where cINdIN and dINdIN indegrees (independent variables are $I}_{cIN>13$ and ${I}_{dIN}$) correlate very strongly with the period of swimming:
where $T$ is the period. The coefficient of determination ${R}^{2}=0.96$.
We used this linear regression model to predict firing period for 200 new connectomes (100 probabilistic, 100 anatomical). The accuracy estimated using the coefficient of determination is ${R}^{2}=0.94$ (Figure 5C). We were therefore able to predict with good accuracy a key characteristic of the network’s behaviour based only on its connectivity. Note that this prediction is universal, since it does not require knowledge of whether the connectome was generated using the probabilistic or anatomical model.
Why is inhibition from cINs stronger relative to excitation from dINs, and therefore swimming slower, in connectomes generated by the probabilistic model? This is a difficult question to answer completely, but much of the difference is due to the fact that anatomical connectomes have more cINs that receive fewer than 13 connections from dINs and are thus inactive during swimming (anatomical model: 168 ± 11 inactive cINs, N = 100; probabilistic model: 101 ± 8 inactive cINs, N = 100). Although the mean dINcIN indegree is very similar between anatomical and probabilistic connectomes (and above the threshold of 13), the variance is much higher in the anatomical case (Figure 5D).
Therefore, in anatomical connectomes there are more inactive cINs. The underlying reason for this difference in variance is that in the anatomical model neurons have randomly chosen dendritic extents, sampled from the distribution of experimentally measured dendrites (see Appendix 1). This means that some neurons have small dendrites and receive very few connections, while others have large dendrites and receive very many connections. In the probabilistic case this detail is lost, as all incoming connections to a neuron are chosen completely independently of each other.
While we can explain the quantitative difference between anatomical and probabilistic models, this difference clearly illustrates that there are potential pitfalls in applying the probabilistic approach to a particular biological question, and it must be used with caution. In this specific case, there is a problem because the reduction in dIN to cIN indegree variance produced by the generalization process used to generate the probabilistic connectomes has asymmetric consequences. The decreased number of cINs failing to fire because of weak excitatory input (low indegree number) is not balanced by the effect of reducing the number of cINs with very strong excitatory input (high indegree number). This is because, above a threshold input strength, cINs only fire a single spike per cycle (see the section entitled "Reliability of cIN spiking depends on their RCcoordinate"); changing the level of excitation above the threshold value does not alter this. The result of this asymmetry is the overall increase in the number of cINs firing with the probabilistic model and hence the lengthened cycle period. To offset this consequence of the probabilistic approach, we therefore reduced the strength of cIN to dIN inhibition (from 0.435 to 0.2 nS). As predicted, this reduced the cycle period to a range overlapping the distribution produced by anatomical connectomes and matches periods seen in the real swimming behaviour (see orange histogram of swimming periods for the modified probabilistic connectomes in Figure 5A).
A core dINcIN subnetwork can generate swimming
The probabilistic approach allows us to test the reliability of network function after removal of selected connections. As an illustration, we considered a subnetwork comprising only the sensory pathway (which is not active during swimming), and dIN and cIN CPG neurons. We excluded aINs and mns simply by setting the probability of connections to and from them to zero. Figure 6 shows one simulation of the functional model containing only this subnetwork. Figure 6A,D shows examples of voltage dynamics for individual dIN and cIN neurons on the right and left body sides, respectively; Figure 6B,C shows raster plots of spiking activity for all neurons on the right and left sides of the body, respectively.
The brown and light blue dots in Figure 6B,C show a typical pattern of antiphase, leftright swimming activity in the dINcIN subnetwork. We found that in 100 independent simulations (with different reduced network connectomes) swimming activity was generated that was similar to that in Figure 6. The swimming period in these simulations was 57 ± 0.9 ms. These values are again within the physiological range observed in experimental recordings of swimming.
Previous experiments have shown that the swimming CPG includes dINs, cINs and aINs (Roberts et al., 2010). However, it is known that aINs have a low probability of firing during swimming, suggesting that their contributions during swimming are minimal and their role in the network is still unclear (Li et al., 2004). Our simulation results confirm these experimental findings by showing that the dINcIN subnetwork generates reliable swimming.
Removal of commissural connections allows rhythmic firing on the stimulated body side
Experiments have revealed that an isolated side of the tadpole spinal cord without commissural connections can generate regular rhythmic spiking activity in motoneurons, with period that is lower than that of swimming (Soffe, 1989).
Once again, the probabilistic model readily allowed us to simulate these experimental findings by setting the probability of commissural connections from cINs and dlc neurons to zero to disconnect the two body sides. This is equivalent to a sagittal midline lesion experimentally. Figure 7A shows a raster plot of steady oscillatory spiking in motoneurons (green) and dINs (brown), demonstrating that the rhythmic activity was maintained and stable.
It is important to note that the mechanism that generates this singlesided rhythm is different to that which generates swimming. In swimming, inhibition from cINs causes contralateral dINs to fire postinhibitory rebound spikes. In the case of separated body sides there is no cIN input to dINs, and the only other inhibitory CPG neurons, the aINs, are inactive. Instead, the rhythmic activity is caused by feedback NMDA excitation within the dIN population, as has been previously observed experimentally (Li et al., 2010) and in modelling (Hull et al., 2016). Within one simulation dINs fell into a number of different groups, based on their spiking period. In most simulations, the majority of dINs spiked rather quickly, with period approximately 24 ms (41 ± 16 dINs, N = 100 connectomes), while most of the remaining dINs spiked with approximately double this period, approximately 53 ms (24 ± 4 dINs). A much smaller group fired twice as slowly again, with a period of approximately 101 ms (2 ± 3 dINs). Figure 7B makes these groups clear, by showing the same set of dIN spikes as Figure 7A but with the neurons sorted according to firing rate. Interestingly, motoneurons tended to fire inphase with the intermediately sized group of dINs that spiked at approximately 53 ms (as shown in Figure 7A), although in some simulations some motoneurons did also spike inphase with the faster group of dINs; further investigation is required to understand why more mns are not able to fire with the dINs in this group.
We have no direct experimental recordings of dINs following separation of the two body sides, only ventral root recordings showing motor neuron activity. In these experiments (Soffe, 1989), it was found that singlesided rhythmic activity was significantly faster than that seen during swimming (initial average cycle period 60 ms vs 43 ms). This was also the case with our simulations, where most mns spiked at approximately 53 ms in the singlesided cases, versus approximately 69 ms in normal swimming. From our results, we predict that recordings from dINs during singlesided rhythm generation would reveal a relatively large group of dINs that spike much more quickly than ventral root activity, and another much smaller group of dINs that fire much more slowly.
Reliability of cIN spiking depends on their RCcoordinate
Experiments have shown that during swimming the reliability of spiking of some neuron types can vary from cell to cell (Soffe, 1993; Li et al., 2007). In simulations of connectomes generated by the anatomical model approximately 50% of cINs fire reliably, whereas in connectomes from the probabilistic model approximately 70% of cINs were reliable. Other cINs were either completely inactive or only fire on some swimming cycles. We investigated the cause of this unreliability by analysing the probabilistic model.
In the functional model, for each pair of cell types, the mean value of synaptic strength was selected in line with experimental data (Roberts et al., 2014) and randomised by addition of the Gaussian random variable with zero mean and relatively small variance (see Appendix 2, Synaptic Currents). In the case of synchronous bombarding, the total input to the neuron depends on both the connection strength and the number of incoming connections, therefore, the degree is an important measure. For the reliability study, we approximate the total input to cIN by the mean dIN to cIN connection strength multiplied by the mean indegree from dINs to cINs, because dIN spike reliably and synchronously during each swimming cycle.
From simulations of 100 different connectomes, we found that the probability that a cIN spikes reliably depends on the dINcIN indegree (${I}_{dIN}$). If ${I}_{dIN}>15$ then a cIN fires once on each swimming cycle, approximately in phase with dINs and mns on the same side; we call this a 'reliable' cIN. If ${13\le I}_{dIN}\le 15$ then firing is irregular, meaning the cIN fires approximately inphase with dINs and mns but on only some swimming cycles; we call this an 'unreliable' cIN. Those cINs that have ${I}_{dIN}<13$ do not fire at all during swimming.
The probabilistic model allowed us to calculate the expected dINcIN indegree as a function of its rostrocaudal position (Figure 8A). Note that this result was based only on analysis of the general probability matrix, not individual connectome realisations. The relationship allowed us to hypothesise: (1) it is likely that rostral cINs will not fire; (2) it is likely that cINs with RCcoordinate near 900 µm are unreliable, and (3) it is likely that caudal cINs will fire reliably.
To confirm these hypotheses in the model we used the results of 100 spiking simulations to calculate the probability that a cIN in a certain position will fire reliably. In Figure 8B we show the reliability proportion (the fraction of simulations where the cIN fires reliably) vs RC coordinate. From this, it was clear that cINs at more rostral positions have a significantly lower probability of reliable spiking than cINs in more caudal positions. Using a linear regression model, we determined a strong correlation between the cIN reliability proportion ($x$) and the average dINcIN indegree (y) given by the linear relationship $y=0.07\bullet x0.4$ (Figure 8C). Note that there is currently not enough experimental data about the reliability of cIN spiking during swimming in vivo to say whether the level of cIN reliability in our simulations was realistic. However, our general results from the probabilistic model suggest that it is important that any experimental measures of cIN spiking reliability (or that of other neuron types) should take into account the rostrocaudal position of the measured cell.
Ascending axons of dINs are important for swimming
It is a defining feature of dINs in the tadpole that they all have a descending axon, but some dINs which are located more rostrally have a second axon growing in the ascending direction (Borisyuk et al., 2014; Roberts et al., 2014). Simplified computational models (Li et al., 2006; Wolf et al., 2009) have shown that the swimming activity fails to selfsustain unless some excitatory interneurons have ascending connections. We used the probabilistic model to further clarify the role of ascending dIN axon branches, taking advantage of the fact that our new model allows us to run large numbers of simulations and to study the generalised connection structure. Using the probabilistic model, we removed all ascending connections from dINs and generated a modified adjacency matrix (connectome), which we then used to simulate spiking activity.
Figure 9A shows the indegrees for the dIN subnetwork (i.e. the number of incoming connections to each dIN from other dINs) for the standard connectome (black) and one lacking ascending dIN axons (red). In the figure, the horizontal and vertical axes show the indegrees and the RCcoordinate of dINs, respectively. We consider here only rostral and midbody dINs in the range of RCcoordinates from 500 to 1400 µm; more caudal dINs do not receive any synapses from ascending dIN axon collaterals, so the indegrees are the same for both connectomes.
From Figure 9A it is clear that the dIN indegrees in both cases are similar in the middle body part but are increasingly different for neurons in the rostral part. For the modified connectome, the indegree (red dots) decays to zero in a linear way as the RCcoordinate approaches 500 µm because the dINs in the most rostral locations have a few if any connections from descending axons. As a result, the most rostral dINs in the modified connectome can only fire due to electrical coupling between dINs, resulting in the appearance of some unusual patterns of spiking activity not observed experimentally. We repeated 100 simulations of the functional model after removing ascending dIN connections. The resulting spiking activity patterns can be divided into three cases:
Case 1 (63/100): In most simulations, the swimming activity was initiated but failed to persist. Swimming failures begin with rostral dINs failing to spike due to reduced excitatory drive from other dINs (ascending dIN connections are missing); this reduced excitation from the rostral dINs prevents slightly more caudal dINs from firing, and so on, as can clearly be seen in Figure 9B. This result is in line with previous modelling that showed that feedback excitation is a mechanism that contributes to generating persistent motor activity in a simpler model (Li et al., 2006).
Case 2 (36/100): In 33 of 36 simulations one side only was active. In 3 of 36 simulations both sides were rhythmically active for the total length of the simulation but they do not fire in antiphase. The pattern of spiking activity on one side is similar to the one shown in Figure 7A.
Case 3 (1/100): Only one simulation generated sustained swimming alternating firing between left and right sides, but the period of the oscillations was shorter than for the standard connectome (50 ms).
Discussion
The study of neuronal connectivity is a challenging problem in contemporary neuroscience. One popular and effective method for finding cortical connectivity involved detailed tracing of a small number of individual neurons of each identified type, and then using estimates of the number of location of the different cell types to estimate complete connectivity (Binzegger et al., 2004). Recent development of new brain imaging techniques allows generation of 3D images of single neurons, tracing their connections and, for example, making progress towards a complete Drosophila connectome (Lin et al., 2015; Shih et al., 2015; National Center for Highperformance Computing and National Tsing Hua University, 2009). Similar progress has been made by combining molecular, anatomical and physiological techniques to find the neuronal cell types, and connections between them, in mouse retina (Seung and Sümbül, 2014; Kim et al., 2014). Computational modelling has been successfully applied to find a sensorimotor connectome in larval Zebrafish (Stobb et al., 2012). In this paper, available neurobiological data have been used to describe neuronal cell types and formulate a stochastic model of connectivity, which was studied using a graph theory approach.
It is known that brain development involves multiple stochastic processes and that, in most species, individuals’ connectomes are different (Seung, 2012). Despite differences in connectivity, most individuals under normal conditions are able to demonstrate similar functionality. This means that different connectomes include sufficient key structural features to produce a common repertoire of functionalities and behaviours. What are the key connectivity properties that define the network functionality?
Motivated by this question, we derive a probabilistic model of connectivity in the Xenopus tadpole CNS (caudal hindbrain and spinal cord) to study the relationship between the structure and function of the network. To derive the probabilistic model we generate 1000 connectomes using a biologically realistic anatomical model based on the ‘developmental’ process of axon growth (Li et al., 2007a; Borisyuk et al., 2011; Borisyuk et al., 2014; Roberts et al., 2014). A similar approach to generating connectivity from a developmental process was used by Bauer et al. (2014); in this case, a reactiondiffusion model was applied to generate connectivity in a network of excitatory and inhibitory neurons with winnertakesall functionality.
Using a universal ordering of neurons in the tadpole, we have calculated the probability of connection from each neuron ( i ) to neuron ( j ) as the frequency at which a connection exists among the thousand generated connectomes. In this way, our probabilistic model ‘generalizes’ structural properties of networks produced by the anatomical model.
Using the probabilistic model, we can generate an adjacency matrix representing a particular realisation of neuronal connectivity. Mapping the adjacency matrix to a functional model of spiking neurons of HodgkinHuxley type enables us to simulate spiking activity. We compare these simulations of the functional model to the experimental results on swimming initiated by skin touch. All generated adjacency matrixes (connectomes) mapped to the functional model generate similar swimming activity. It seems, then, that the probabilistic model contains some fundamental features of the network connectivity (‘proper structure’) which ensure correct functioning of the system. For example, experimental recordings show that apparentlypathological activity (synchrony) can sometimes appear soon after swimming initiation: the two body sides spike synchronously during several cycles before then switching to normal antiphase swimming activity (Li et al., 2014). This synchronous activity appears also in model simulations with connectivity generated by both the anatomical and the probabilistic models. However, the number of synchronously firing neurons is significantly reduced in probabilistic connectomes.
A second type of apparently pathological activity is the additional firing of some dINs near the middle of the swimming cycle (midcycle dINs) (Li et al., 2014). Midcycle dINs appear in model simulations with both anatomical and probabilistic connectivity. However, the number of such midcycle dINs is significantly reduced in probabilistic connectomes: 0.8 and 6.3 for probabilistic and anatomical connectomes, respectively (average according to swimming cycles and 100 simulations).. These results suggest that synchrony and midcycle dINs arise from connectivity imperfections and that the generalised connectivity encapsulated in the probabilistic model improves on the imperfection of some individual realisations.
To design the probabilistic model, we use a minimalistic approach. We use the assumption that directed connections are represented by the matrix of independent Bernoulli random variables. One of the strengths of this approach is that it allowed us to analytically calculate some of the graph’s characteristics (the mean and standard deviation of in and outdegrees, heterogeneity coefficients) directly from the probability matrix, without considerations of a particular (generated) connectome. In the case of the anatomical model, we can only compute graph characteristics for a connectome realization. Here, we study how these characteristics relate to particular functional properties of the network. For instance, the average in and outdegrees were used to predict the swimming period and to find the positions of reliably firing cINs.
The assumption that the Bernoulli variables are independent is a significant limitation of the probabilistic model. One way to overcome this limitation might be the use of more sophisticated probabilistic processes where the random variables corresponding to different connections become dependent (e.g. random Markov field approach).
Computational modelling of the tadpole spinal cord reveals the fundamental features of neuronal connectivity that are responsible for robust swimming generation. Unlike simpler organisms such as C. elegans, tadpoles have the potential for significant variation between individuals in terms of neuronal connectivity, as a result of the large number of random processes involved in their development. Despite this variation, the behaviour of individuals is approximately the same, suggesting some fundamental organisational principles common across the species. We adopt the philosophy that, for tadpoles at least, there is a theoretical ‘perfect’ version of the nervous system with individual random variations from this ideal. Although, the probabilistic model arises from ‘averaging’ of many anatomical connectomes, this model still generates connectomes that reliably swim and this fact presumably reflects the fundamental organisational principles of the system. An interesting property of connectomes generated by the probabilistic model is that their anatomical and functional characteristics are considerably less variable than those generated by our anatomical model (and on whose properties the probabilistic model was based). We hypothesise that due to the ‘averaging’ process of the probabilistic model, the connectomes generated from it are closer to the theoretical ‘ideal’ network. Some characteristic features of the connectivity are not clear from an individual realisation, but become evident from the probabilistic model. For example, the shape of degree distributions as a function of cell position cannot be clearly seen from analysing an individual connectome – these shapes are irregular. They are much clearer when calculated directly from the probabilistic model itself. In addition to this, connectomes generated by the probabilistic model generate spiking activity that is considerably less variable and ‘messy’ than anatomical connectomes, which makes it easier to see and quantify phenomena such as irregularly spiking cINs.
Finding neural connection probabilities under biological constraints is a difficult problem. In the case of the tadpole spinal cord, the system is simple enough that it is possible to reconstruct biologically realistic connectivity (Roberts et al., 2014) (an anatomical connectome) and to define neuronal connection probabilities (probabilistic model). We believe that this is a promising general approach that could be used beyond the particular case of tadpoles. Similar probabilistic approaches have been used for modelling the development of neural networks using limited experimental data (Binzegger et al., 2004; Zubler and Douglas, 2009). Another possible approach for finding connection probabilities is to minimize an appropriate cost function which reflects both anatomical and functional properties. A combination of these approaches has been used in pilot studies that aim to incorporate into the tadpole connectome a new group of neurons recently found in the hindbrain (Buhl et al., 2015). We will report this result in a separate publication.
Conclusion
We study the structure and function of the spinal cord neuronal network using experimental data and computational modelling. Our anatomical model generates multiple highly variable and nonhomogeneous connectomes and to deal with this large and complex data we design a very simple mathematical metamodel expecting that this new probabilistic model will reflect (generalise) structural properties of anatomical connectomes and show proper functioning.
The crucial question is: ‘Can probabilistic connectomes produce swimming’? The answer to this is not obvious. Our earlier paper (Li et al., 2007a) showed that a graph of connections based on probabilities derived from small number of pairwise recordings provides swimming in about 60% of cases only. On the other hand, this new study shows that probabilistic connectomes that include some of the structure of anatomical connectomes reliably swim in all cases. Thus, we can derive an important conclusion that the two properties of the probabilistic model inherited from anatomical connectomes: (1) position of neurons along the rostrocaudal coordinates and (2) frequency of connection appearance, are sufficient for swimming generation.
Also, it is easy to use the probabilistic approach to generate connectomes compared to the need to ‘grow’ them using the anatomical model: all traditional characteristics of the connectivity graph can be calculated directly from the probability matrix without consideration of particular connectomes. Some characteristics of the probabilistic connectomes (e.g. the mean of in and outdegrees) coincide with equivalent characteristics of the anatomical connectomes but some differ (e.g. the variances of in and outdegrees are significantly smaller for probabilistic connectomes). Although there are some differences between the behaviour of anatomical and probabilistic connectomes, even studying these differences can provide important insights into the relationship between the structure and function of the network. Our investigation in the reasons underlying a difference in swimming frequency between the two types of connectome (see result section) is an example of this, where we found that it was the degree of variance of cIN indegree from dINs that largely caused the difference. It would have been hard to observe this interesting phenomenon without having the probabilistic model (where indegree variance is much lower) to compare with the anatomical one.
The probabilistic model provides a different way to look at the information generated by the anatomical model. It is grounded in the previous anatomical model as the anatomical model is grounded in the biological anatomy. It provides a different perspective on data generated by many anatomical models, and it is this different perspective that makes the probabilistic model an advance.
Appendix 1
Anatomical model
The anatomical model generates the complete neural connectivity in the spinal cord using a ‘developmental’ approach that mimics axon growth (Borisyuk et al., 2014). The axon growth model is based on numerous anatomical data of Xenopus tadpole spinal cord neurons. Our attempt is to include in the model as much biological information as possible. Here, we give a brief description of the model, but full details can be found in (Borisyuk et al., 2014; 5, 4).
The model spatial structure
Tadpole’s spinal cord is approximated as a 2D rectangular plate where neuronal bodies, dendrites and axons are located (Figure 1 in the main text). The third dimension (i.e. the thickness of the spinal cord ‘tube’) was ignored as this is very thin (10 µm thick). In the model description, variables $x$ and $y$ correspond to the rostrocaudal distance (RC) from the midbrainhindbrain border and the dorsoventral distance (DV) from the ventral midline, respectively. Positive (Negative) values of $y$ correspond to positions on the left (right) side of the body. We consider a limited area of the spinal cord, where $\left(x,y\right)\in \left[500\text{}\mu \mathrm{m},2000\text{}\mu \mathrm{m}\right]\times \left[145\mu \mathrm{m},145\mu \mathrm{m}\right]$.
Axon growth model
We describe the axon growth using discrete time iteration map (with time step 1 ms). The map is described by three variables $({x}_{n},{y}_{n},{\theta}_{n})$, where ${x}_{n}$ represents RC coordinate, ${y}_{n}$ the DV coordinate and ${\theta}_{n}$ the growth angle of the axon at each time step $n(n=\mathrm{0,1},\dots ,N)$. The map for the growth angle depends on a ‘stiffness’ term, which is the tendency of the growth angle to grow straight, and by the influence of environmental cues (according to chemical gradients), which deviate the growth cone from a straight path. The chemical gradients functions ${G}_{RC}(x,y)$ and ${G}_{DV}(x,y)$ depend on the current position of the axon and that will determine the change of the growth angle at each time step on the RC and DV axis, respectively. Additionally, a uniform random variable ${\u03f5}_{n}$ is included to provide an additional degree of freedom at each time step. The map is described by the following equations:
Here, the elongation parameter$\mathrm{\Delta}=1\mu m$; ${\u03f5}_{n}$ uniformly distributed in the interval $\left[\alpha ,\alpha \right]$. To start the axon growth simulation we assume that the axon initial positions $({x}_{0},{y}_{0})$ coincide with the soma positions. The initial value for the growth angle ${\theta}_{0}$ and the axon length $L=N\bullet \mathrm{\Delta}$ (which will determine the number of iterations of the map) are randomly selected from the distribution of experimentally measured initial angles and axon lengths.
The gradient functions ${G}_{RC}$ and ${G}_{DV}$ depend on various parameters that describe the properties of the chemical gradients in the 2D space. Previous studies revealed that the axons of the same neuron type tend to grow in specific regions of the spinal cord, suggesting that such axons could be controlled by the same gradients. Thus, the parameters of the gradient functions in the model were selected to reproduce the statistical properties of the axons of each specific cell type separately. Some of these parameters were selected according to general biological knowledge on the distribution and properties of the chemical gradients. The remaining parameters were estimated using an optimization technique that minimizes a custom cost function that measures the similarity between statistical properties of simulated and experimentally measured axons. For a detailed description of the gradient functions ${G}_{RC}$ and ${G}_{DV}$, the cost function and the optimization technique see (Borisyuk et al., 2014.
Anatomical and physiological studies in tadpoles have revealed that there are physical constraints to the axons growth, and these have been added in the anatomical model. In particular, longitudinal barriers delimited by DV coordinates $y=\pm 125$ delimit the area where axons of all neuron types except RBs grow (called marginal zone). Coordinates $y\in \left(\mathrm{127,137}\right)$ on the left body side ($y\in \left(137,127\right)$ for the right body side, respectively) delimit the region where RB axons can grow (called dorsal tract area).
Branching and commissural axons
The axon of tadpole spinal neurons typically splits into two branches during its growth. We call primary axon the branch that starts from the soma and continues to grow following the growth direction before the branching point. Secondary axons grow from the branching and grow towards the opposite direction of the primary axon (rostrocaudally).
The RC direction of primary axon of each neuron type is known from anatomy of neurons. In the model we thus consider growth of primary axons and secondary axons as two consecutive processes. All primary axons grow first, and after that secondary axon growth starts at a specified branching point. The coordinates of branching points are selected randomly from the distribution of available experimental data for each neuron type.
The axon of commissural neurons (dlcs and cINs) starts to grow on the cell body side and then rapidly navigates on the opposite body side after crossing the floor plate due to the influence of strong DV gradients. After crossing, DV gradients change their sensitivity and become weak, and the axon starts to deviate towards ascending direction. Secondary axons are positioned on the contralateral side and grow in the descending direction. At the beginning, commissural neurons grow in the ventral direction according to the axon growth equations with specially adjusted parameter values. After crossing the boundary of the ventral plate on the opposite side the axon growth is described by the same equations but with another regular set of parameter values (for details see (Borisyuk et al., 2014)).
Synapses and full connectome
Dendrites are assumed to be fixed bars extending dorsoventrally. Thus, each dendrite is represented by a pair of DV positions, one corresponding to its lower (ventral) bound and on to its upper (dorsal) bound. For each neuron, cell body positions and dendritic bounds are sampled from the distributions of experimentally measured data for each specific cell type.
In model simulations the number of neurons of each type is fixed and it is the same for both sides of the body (see methods section in the main text). In reality, although biological data are limited, total numbers of spinal neurons and the population sizes of individual neuron types do not appear to vary greatly between animals (perhaps ±10–15% at most) at this early stage of development. Any variation that there is in numbers is small compared to the differences from the previous developmental stage and the following stage, both of which swim in the same way.
For each pair of pre and post synaptic neurons, a connection is generated whenever the presynaptic axon crosses the dendritic bar of the postsynaptic neuron with some probability. Such probability depends on pre and post synaptic cell types and was estimated from various experimental pairwise recordings (for details see (Borisyuk et al., 2014)). During the initial precrossing stage we assume that the axon of commissural neurons cannot produce synapses.
Since the model uses several number of random variables, each simulation of the anatomical model generates a different connectome.
Appendix 2
Functional model
The probabilistic and anatomical connectomes provides detailed information on connectivity in the spinal cord – specifically a list of synaptic connections between neurons. We use this information to build a functional model that simulates the spiking activity of spinal cord neurons. This allows us to study one of the fundamental problems of neuroscience: the relationship between connection structure and functionality.
Overview
To simulate the activity in the generated connectomes we represent each cell as a single compartment conductance based neuron of HodgkinHuxley type. The equation governing the membrane potential (V) for neuron i is:
The capacitance $C$ of all neurons is $10pF$, which corresponds to a density of $1.0\mu F/c{m}^{2}$ for a total surface area of ${10}^{5}c{m}^{2}$. The terms ${I}_{lk}$, ${I}_{Na}$, ${I}_{Kf}$,${I}_{Ks}$ and ${I}_{Ca}$ represent transmembrane currents mediated by different ions, respectively: nonspecific leak, sodium, fast potassium, slow potassium and calcium. The terms ${I}_{syn}$ and ${I}_{gj}$ represent the summed inputs from chemical synapses (${I}_{syn}$) and gap junctions (${I}_{gj}$), while ${I}_{ext}$ is an externallyinjected current. Although the different neuron types in the tadpole spinal cord have different electrophysiology, for simplicity we use the model of a motoneuron from for most model neurons, as this shows characteristics (e.g. repetitive firing in response to injected current) that are broadly shared by all of the neuron types. The exception to this is dINs, which have special properties such as only firing a single spike in response to current injection and the ability to fire postinhibitory ‘rebound’ spikes. Model dINs differ from nondINs in the following ways:
The parameter values governing the membrane properties are different (Appendix 2—table 2).
Only dINs contain a calcium current. For nondINs, we set ${I}_{Ca}=0$.
Only dINs make gap junction connections (and only with other dINs).
Membrane channels
The leak, sodium and potassium channel currents are given by the following equations:
The parameters ${E}_{lk}$, ${E}_{Na}$ and ${E}_{K}$ give the reversal potential for the leak, sodium and potassium channels respectively, and the parameters ${\stackrel{}{g}}_{lk}$, ${\stackrel{}{g}}_{Na}$, ${\stackrel{}{g}}_{Kf}$and ${\stackrel{}{g}}_{Ks}$ give their maximum conductances (these parameter values are given in Appendix 2—table 1). The gating variables $h$, $m$, ${n}_{f}$ and ${n}_{s}$ are governed by equations (Equation 6; Equation 7; Equation 8; Equation 9), where $X=h,m,{n}_{f},{n}_{s}$.
The values of the parameters A, B, C, D and E in the functions ${\alpha}_{X}\left(V\right),{\beta}_{X}\left(V\right)$ were taken from Sautois et al. (2007) for nondINs and from (Roberts et al., 2014) for dINs, and are shown in Appendix 2—table 2.
As in (), model dINs contain a calciummediated current which is modelled according to the GoldmanHodgkinKatz equation. This current is calculated as:
Here, ${p}_{ca}$ is the permeability of the membrane to calcium ions (analogous to maximum conductance) and $z$ is their ionic valence (+2). ${S}_{in}$ and ${S}_{out}$ are the concentration of calcium in and outside of the cell, respectively. F is Faraday’s constant, and R is the ideal gas constant, while T is the temperature in Kelvin. Parameters of the calcium current are ${p}_{ca}=14.25c{m}^{3}/ms$, $F=96485C/mol$, $R=8.314J/\left(Kmol\right)$, $T=300K$, ${\left[C{a}^{2+}\right]}_{i}={10}^{7}mol/c{m}^{3}$, ${\left[C{a}^{2+}\right]}_{o}={10}^{5}mol/c{m}^{3}$. Finally, ${h}_{Ca}$ is the gating variable associated with the calcium current, which is governed by the standard gating equations (Equation 6; Equation 7; Equation 8; Equation 9) – although note from Appendix 2—table 2 that two different sets of parameters are used for this equation based on whether the membrane potential is above or below −25 mV.
Synaptic currents
The synaptic current that arises in a neuron is the combination of three different subtypes of synaptic receptor: excitatory AMPA and NMDA and inhibitory glycine:
Each synaptic current is calculated using the following equation, where $X=ampa,nmda,inh$:
Here, ${\stackrel{}{g}}_{i,j}^{X}$ is the maximum conductance (‘strength’) of synaptic connection of type X from neuron $i$ to neuron $j$. If the connectome does not include a connection from $i$ to $j$ then ${\stackrel{}{g}}_{i,j}^{X}=0$, otherwise it is selected according to the type of the pre and post synaptic neurons, based on paired recordings. Presynaptic neuronal type determines the synapse type $X$. Inhibitory neuron types are cIN and aIN; excitatory ones are the remaining cell types. The synaptic strengths used in the model are typically the ones given in Appendix 2—table 3, except for few values that were modified to match the physiology of neurons and synapses (details are given in (Roberts et al., 2014)). Specifically, the maximal conductance of AMPA synapses from RBs to dli neurons are set to $8nS$, the maximal conductance of AMPA synapses from dINs to aINs are set to $0.1nS$, the maximal conductance of NMDA synapses from dIN to dINs are set to $0.15nS$, and the maximal conductance of NMDA synapses from RB to dlc are set to 1nS..
The set ${S}_{j}\left(t\right)$ contains the times of all the spikes that neuron $j$ has fired up to the current time $t$. Each spike generates a postsynaptic current (PSP) that rises according to the time constant ${\tau}_{o}^{X}$ and decays according to ${\tau}_{c}^{X}$. The normalizing constant ${\mathrm{\Delta}}_{X}$ is set such that the peak of the sum of the exponentials is 1, meaning that following a spike the conductance rises to a maximum of ${\stackrel{}{g}}_{i,j}^{X}$. The values selected for the time and normalizing constants are given in Appendix 2—table 3 and they are based on previous modelling (Roberts et al., 2014. To mimic synaptic strength variability, Gaussian noise with standard deviation 5% of the mean was added to the maximum conductance of each individual synapse.
The synaptic delay between two neurons, ${\delta}_{i,j}$, consists of a constant and distancedependent part:
Here, ${P}_{i}$ and ${P}_{j}$ are the positions of neurons $i$ and $j$ along the rostrocaudal axis, ${\delta}_{C}$ is the constant delay and ${\delta}_{D}$ is the speed of synaptic transmission. We set ${\delta}_{C}=1ms$ and ${\delta}_{D}=0.0035ms/\mu m$.
Finally, the function ${f}_{X}\left(V\right)$ determines how the synaptic current depends on the postsynaptic voltage. For AMPA and inhibitory synapses this has a simple linear (Ohmic) form:
Where $X=ampa,inh$ and ${E}_{X}$ is the equilibrium (reversal) potential of the synapse type (Appendix 2—table 3). As a result of magnesium block, NMDA synapses have an additional nonlinear voltage dependence, which we include by adding a sigmoidal scaling term to ${f}_{NMDA}$:
Gap junctions
Descending interneurons (dINs) are electrically coupled to other nearby dINs via gap junctions. For dINs only, the gap junction current is calculated using a simple Ohmic relationship:
Here ${G}_{i}$ is the set of indexes of all dINs that are on the same side of the body as neuron $i$ and are located within ${D}_{gj}$ of neuron $i$ on the rostrocaudal axis, where we set ${D}_{gj}=100\mu m$. The parameter ${\stackrel{}{g}}_{gj}$ gives the conductance of gap junctions, and we use the value ${\stackrel{}{g}}_{gj}=0.2nS$.
References
 1
 2

3
A quantitative map of the circuit of cat primary visual cortexJournal of Neuroscience 24:8441–8453.https://doi.org/10.1523/JNEUROSCI.140004.2004

4
Modeling the connectome of a simple spinal cordFrontiers in Neuroinformatics 5:20.https://doi.org/10.3389/fninf.2011.00020
 5

6
Sensory initiation of a coordinated motor response: synaptic excitation underlying simple decisionmakingThe Journal of Physiology 593:4423–4437.https://doi.org/10.1113/JP270792

7
Complex brain networks: graph theoretical analysis of structural and functional systemsNature Reviews Neuroscience 10:186–198.https://doi.org/10.1038/nrn2575
 8
 9

10
Flycircuit DatabaseAccessed April 20, 2018.

11
Unified index to quantifying heterogeneity of complex networksPhysica A: Statistical Mechanics and its Applications 387:3769–3780.https://doi.org/10.1016/j.physa.2008.01.113
 12
 13

14
The brainstem reticular formation is a smallworld, not scalefree, networkProceedings of the Royal Society B: Biological Sciences 273:503–511.https://doi.org/10.1098/rspb.2005.3354
 15
 16
 17

18
Primitive roles for inhibitory interneurons in developing frog spinal cordJournal of Neuroscience 24:5840–5848.https://doi.org/10.1523/JNEUROSCI.163304.2004
 19
 20
 21
 22

23
Persistent responses to brief stimuli: feedback excitation among brainstem neuronsJournal of Neuroscience 26:4026–4035.https://doi.org/10.1523/JNEUROSCI.472705.2006

24
Automated in situ brain imaging for mapping the Drosophila connectomeJournal of Neurogenetics 29:157–168.https://doi.org/10.3109/01677063.2015.1078801

25
Principles of rhythmic motor pattern generationPhysiological Reviews 76:687–717.https://doi.org/10.1152/physrev.1996.76.3.687
 26

27
How neurons generate behavior in a hatchling amphibian tadpole: an outlineFrontiers in Behavioral Neuroscience 4:16.https://doi.org/10.3389/fnbeh.2010.00016
 28

29
Role of typespecific neuron properties in a spinal cord motor networkJournal of Computational Neuroscience 23:59–77.https://doi.org/10.1007/s1082700600191
 30
 31

32
Connectomicsbased analysis of information flow in the Drosophila brainCurrent Biology 25:1249–1258.https://doi.org/10.1016/j.cub.2015.03.021
 33
 34
 35

36
The human connectome: A structural description of the human brainPLoS Computational Biology 1:e42.https://doi.org/10.1371/journal.pcbi.0010042

37
The small world of the cerebral cortexNeuroinformatics 2:145–162.https://doi.org/10.1385/NI:2:2:145
 38
 39

40
The rich club of the C. elegans neuronal connectomeJournal of Neuroscience 33:6380–6387.https://doi.org/10.1523/JNEUROSCI.378412.2013

41
Neural development features: spatiotemporal development of the Caenorhabditis elegans neuronal networkPLoS Computational Biology 7:e1001044.https://doi.org/10.1371/journal.pcbi.1001044

42
Structural properties of the Caenorhabditis elegans neuronal networkPLoS Computational Biology 7:e1001066.https://doi.org/10.1371/journal.pcbi.1001066
 43

44
A framework for modeling the growth and development of neurons and networksFrontiers in Computational Neuroscience 3:25.https://doi.org/10.3389/neuro.10.025.2009
Decision letter

Ronald L CalabreseReviewing Editor; Emory University, United States
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Structural and functional properties of a probabilistic model of neuronal connectivity in a simple locomotor network" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom, Ronald Calabrese is a member of our Board of Reviewing Editors and the evaluation has been overseen by Timothy Behrens as the Senior Editor.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
This manuscript reports on a modeling approach for building network models based on connectomics data. Using a developmental (anatomical) model of the Xenopus hatchling tadpole spinal cord connectome, the authors generated 1,000 conectomes from which they derive a complete probability matrix of connections in the spinal cord. From this probability matrix (probabilistic model) they then generate ensembles of 100 connectomes (model instances) to simulate based on previous functional models of swimming activity. They find that these instances of the probabilistic model are indeed functional and more robust that models based on developmental model connectomes, which they attribute to lower variance among connectomes for the probabilistic model. They then analyze these model instances to determine functional properties of the swim network, for example showing that the network does not appear to have hub neurons as an organizing principle and explicating how activity in a heminetwork arises. The technique may applicable to connectomics data and provide valuable insights that could benefit both modelers and experimentalists in the field.
Essential revisions:
The reviewers had partially overlapping concerns, which we have tried to combine as much as possible, but some overlap will still be noted in the required revisions below that represent a strong consensus. All minor concerns from the expert reviewers should also be addressed.
1) Both introduction and discussion mention the probabilistic model extracts fundamental structural features. However, the model consists of a celllevel probability matrix, which suggests these structural features may not be so fundamental but specific to the particular assumptions in this circuit model. Furthermore, since the probabilistic model is derived from another model (the anatomical) and not directly from experimental data, it is not clear whether it might just be capturing properties of this anatomical model that don't really generalize to the real system. Moreover, because the probabilistic description is cellwise, it may just reflect an overfitting to the anatomical model. The authors should provide further justification of (1) how well the anatomical model captures the real system, and (2) why/how the probability model is capturing fundamental structure features, despite being at a celllevel granularity. This will also help determine how significant/valid the analyses and predictions made by the model are.
2) The analysis of generative models of neuronal networks is a promising avenue to explore links between development, structure and function. However, we find it difficult to understand the results of this paper relative to the Roberts et al., 2014 work on which it builds. We encourage the authors to either restructure the paper in a way that clarifies not only the differences from the previous model, but the advances in the current work and whether they arise from new computational experiments or the new modeling approach.
Once this clarification is achieved then it will be important to show the significance of some of the conclusions. For example, the fact that all connectomes achieve swimming or "this pathological midcycle spiking is significantly reduced in probabilistic connectomes" need to be contextualized. Since the probability matrix was directly extracted by averaging over anatomical model instances (all of which showed swimming and no pathologies), perhaps it is not surprising that they show similar properties.
3) While the Bernoulli trial formulation in principle admits more analytical analysis, it is primarily used to compute in and outdegrees that could just as easily be computed from direct consideration of the results of the anatomical model. As a result, the probabilistic model is a particularly complex null model that captures many aspects of dynamics and not others, but in ways that are only weakly controlled.
In particular, it is not clear why it is better to analyze the ensemble average properties of individual nodes from the probabilistic model rather than analyze the distribution of results from the previous anatomical model. The strongest observation is in subsection “Functional properties of the model: reliable swimming”, effectively relating the swim period to the variance in indegree of cINs. In the probabilistic model, the use of a binomial distribution implicitly links mean and variance in a way that clearly isn't true of the anatomical model and is shown meaningful impact on the resulting dynamics. This result could be greatly strengthened, if the authors used a random process where variance could be defined separately from the mean and thus this relationship explored explicitly.
4) The last paragraph in the Discussion section makes it seem as if the authors are proposing a novel methodology. However, this novelty is not clear since employing probability matrices to define models is a standard practice: it’s a feature included in several major neural simulators (NEST, Brian, Moose,.…), and many existing models (e.g. on ModelDB) employ prob matrices, either at a populationlevel or celllevel (for small networks and when enough exp data is available).
5) Authors should include more details about the functional model employed since this is crucial for understanding the results (e.g. low connection weights can drastically alter the effect of high connection probabilities). A table with the main parameters used in the functional model should be included. Additionally, authors must share all the model and analysis code  this is common practice nowadays in computational neuroscience to ensure replicability and reproducibility.
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "Structural and functional properties of a probabilistic model of neuronal connectivity in a simple locomotor network" for further consideration at eLife. Your revised article has been favorably evaluated by Timothy Behrens (Senior editor), a Reviewing editor, and two reviewers.
The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:
1) The explanation of why the probabilistic model is an advance over the anatomical model is not yet convincing. Looking at both the response to the reviewers and the argument in the Introduction, it still seems that the same results could have been obtained by consideration of the distribution of networks generated by the anatomical model. The authors must clarify why they emphasize the need to construct the probabilistic matrix as a new step in the analysis process, otherwise averaging the outcome of an existing model and putting into a simple probabilistic framework is a modest and intuitive extension of the previous work without added value.
2) The authors have not fully addressed the issue of why the celllevel probabilistic model is capturing general/fundamental structural features of the tadpole network and not overfitting.
Without addressing these two issues more thoroughly, the impact of the paper rides on its scientific results, connecting the structure of the model to key aspects of its functional output and the authors' emphasis on method is misplaced. Some clarifying points from the expert reviews are presented below.
Reviewer #3:
I don't think the authors have addressed the issue of why the celllevel probabilistic model is capturing general/fundamental structural features of the tadpole network and not overfitting. They argue that it is the result of averaging over 1000 anatomical model connectomes, but these are still a model, e.g. they have a fixed number of neurons for each population, which is not the case for real tadpole networks, and so describing the structure at the celllevel prevents it from being applicable to specimens with a different number of neurons. They also argue that all connectomes achieve reliable swimming, but I fail to see how this is relevant. I think the paper still has value despite this limitation, but I think it should be clearly stated in the paper.
https://doi.org/10.7554/eLife.33281.018Author response
Essential revisions:
The reviewers had partially overlapping concerns, which we have tried to combine as much as possible, but some overlap will still be noted in the required revisions below that represent a strong consensus. All minor concerns from the expert reviewers should also be addressed.
1) Both introduction and discussion mention the probabilistic model extracts fundamental structural features. However, the model consists of a celllevel probability matrix, which suggests these structural features may not be so fundamental but specific to the particular assumptions in this circuit model. Furthermore, since the probabilistic model is derived from another model (the anatomical) and not directly from experimental data, it is not clear whether it might just be capturing properties of this anatomical model that don't really generalize to the real system. Moreover, because the probabilistic description is cellwise, it may just reflect an overfitting to the anatomical model. The authors should provide further justification of (1) how well the anatomical model captures the real system, and (2) why/how the probability model is capturing fundamental structure features, despite being at a celllevel granularity. This will also help determine how significant/valid the analyses and predictions made by the model are.
We agree that this comment is very crucial for understanding the key message of the manuscript and we would like to thank the reviewers for important and fruitful suggestions. Our response and summary of revision related to these comments are explained below. We provide additional description on how the anatomical model captures the real system in the Appendix.
1) The anatomical model captures several important structural properties of the real biological object.
Firstly, in the Conclusion section we have shown that the anatomical model is biologically realistic. The model mimics the growth of axons inside a region of the spinal cord, and parameters of the model are optimized to produce the same statistical characteristics as the experimentally measured axons. Moreover, the model uses numerous biological data (including location of neuron bodies, dendritic bounds, and axon properties) and general anatomical information known from experiments.
Secondly, in the Discussion section we have shown that the connection architecture generated by the anatomical model can be combined with a functional model of spiking units to demonstrate very reliable swimming activity, similar to the one observed experimentally.
Relevant changes in the Introduction (also look at Appendix):
“This anatomical model mimics the realistic growth of axons in the spinal cord. Following biological realism, the axon growth is guided by the concentration of chemical gradients in the spinal cord. The properties of such gradients are controlled by model parameters that have been optimized to produce the same statistical characteristics as real measurements. Other model specifications (including soma positions and dendritic extents) are assigned from the distributions of experimental data and from general biological knowledge”
2) The probabilistic model (matrix) for each pair of neurons gives the probability of directed connection. Each probability is the result of averaging 1000 connectomes generated by the anatomical model. Thus, despite having celllevel granularity, the probabilistic model reflects the biological structure of the real system, because it is generated from the anatomical connectomes.
The probabilistic model generates connectomes and the fundamental feature of these probabilistic connectomes is that the functional model with any probabilistic connectome demonstrates very reliable swimming. Moreover, comparison of swimming firing patterns produced by anatomical and probabilistic connectomes shows that in the case of the probabilistic connectome the firing pattern is closer to the “ideal” swimming pattern. For example, there is a reduction in pathological activities, such as less midcycle dIN firings.
Relevant changes in the Introduction:
“Being derived from multiple biologically realistic (anatomical) connectomes, the probabilistic model reflects the anatomical structure of the biological system.”
“Multiple functional simulations of probabilistic connectomes demonstrated a reliable pattern of rhythmic activity, qualitatively like tadpole swimming and as seen in previous modelling (Roberts.et al., 2014). Thus, the generalised probabilistic model shares structural and functional properties with the real biological object.”
Relevant changes in the Discussion section:
“This synchronous activity appears also in model simulations with connectivity generated by both the anatomical and the probabilistic models. […].These results suggest that synchrony and midcycle dINs arise from connectivity imperfections and that the generalised connectivity encapsulated in the probabilistic model improves on the imperfection of some individual realisations.”
3) We find that in the probabilistic model the variability of in and out degrees is less than in the anatomical model and analyse the consequences of this (comparison of swimming periods in case of probabilistic and anatomical connectomes).
Relevant changes the Introduction:
“Specifically, we found the number of incoming connections (indegree) or outgoing connections (outdegree) of each neuron is higher in anatomical rather than probabilistic connectomes. As a result of this finding, we observed that the period of the rhythm was longer in probabilistic connectomes.”
We appreciate the reviewers’ concern about overfitting. However, the probabilistic model estimates probability of each connection separately and independently from other connection probabilities. Each estimate is based on 1000 cases/observations from the anatomical model.
Relevant changes the Introduction:
“It is important to note that the anatomical model provides a way of generating many different connectomes, such that the random variation observed between generated connectomes has the same statistical properties as measurements taken from different individual animals. This helps to avoid the problem of overfitting the data, as we were able to show that all generated connectomes were able to produce reliable swimming.”
2) The analysis of generative models of neuronal networks is a promising avenue to explore links between development, structure and function. However, we find it difficult to understand the results of this paper relative to the Roberts et al., 2014 work on which it builds. We encourage the authors to either restructure the paper in a way that clarifies not only the differences from the previous model, but the advances in the current work and whether they arise from new computational experiments or the new modeling approach.
Once this clarification is achieved then it will be important to show the significance of some of the conclusions. For example, the fact that all connectomes achieve swimming or "this pathological midcycle spiking is significantly reduced in probabilistic connectomes" need to be contextualized. Since the probability matrix was directly extracted by averaging over anatomical model instances (all of which showed swimming and no pathologies), perhaps it is not surprising that they show similar properties.
This paper describes a new probabilistic model, which is based on our previous work. All our results follow from new modelling and numerous computational experiments. We revised and added the text below to clarify this statement and highlight the significance of our conclusions.
1) In Roberts et al., 2014 we reported both the anatomical model (to generate connectomes) and the functional model (to produce spiking activity of the spinal cord). In the current manuscript, we describe a new probabilistic model, which also generates connectomes. First, we show that the probabilistic model generates connectomes with properties that match the biological object. Second, we demonstrate several important advantages of the probabilistic model in comparison with the anatomical one. Third, we show that the probabilistic model allows us to explain some experimental data. For example, the probabilistic model reveals how the position of cIN neurons is correlated to their firing reliability.
Relevant changes in the introduction:
“Despite the differences between anatomical and probabilistic models, we demonstrate several important advantages of using the probabilistic model in comparison to the anatomical one. For example, we could predict the position of commissural interneurons (cINs) that are active during swimming, which cannot be explained by the anatomical model. Specifically, our simulations show that cINs in rostral positions are less likely to fire reliably than those in caudal positions. Moreover, the probabilistic model allowed us to easily design new computational experiments that helped to clarify the following experimental findings.”
2) Although both the probabilistic and anatomical connectomes generate reliable swimming, the anatomical connectome in many cases produces midcycle dIN spiking and synchrony at the start of a swimming episode which are pathological activities.Therefore, in this respect the probabilistic model is superior to the anatomical one and generates an activity pattern which is closer to the ideal swimming pattern. To better demonstrate this conclusion and provide additional context we have added some results on comparison of motor neuron spike reliability.
Relevant changes in the Discussion section:
“This synchronous activity appears also in model simulations with connectivity generated by both the anatomical and the probabilistic models. […]. These results suggest that synchrony and midcycle dINs arise from connectivity imperfections and that the generalised connectivity encapsulated in the probabilistic model improves on the imperfection of some individual realisations,”
3) While the Bernoulli trial formulation in principle admits more analytical analysis, it is primarily used to compute in and outdegrees that could just as easily be computed from direct consideration of the results of the anatomical model. As a result, the probabilistic model is a particularly complex null model that captures many aspects of dynamics and not others, but in ways that are only weakly controlled.
In particular, it is not clear why it is better to analyze the ensemble average properties of individual nodes from the probabilistic model rather than analyze the distribution of results from the previous anatomical model. The strongest observation is in section 3.2, effectively relating the swim period to the variance in indegree of cINs. In the probabilistic model, the use of a binomial distribution implicitly links mean and variance in a way that clearly isn't true of the anatomical model and is shown meaningful impact on the resulting dynamics. This result could be greatly strengthened, if the authors used a random process where variance could be defined separately from the mean and thus this relationship explored explicitly.
In the current manuscript, designing the probabilistic model, we use the minimalistic approach to construct as simple model as possible. Therefore, we use the assumption that directed connections are represented by a matrix of independent Bernoulli variables. This theoretical assumption allows us analytically to calculate some of the graph’s characteristics (the mean and the variance of in and outdegrees, heterogeneity coefficient) without consideration of generated connectomes. For example, the distribution of indegrees is described by the Poisson binomial distribution and it allows us to use formulas for the mathematical expectation and the variance of indegrees. We have modified Figure. 3A to show the standard deviation of indegree of the probabilistic model.
In case of anatomical model, we can only compute the graph’s characteristics for a particular realization of connections.
We agree that the result can be strengthened in case of using the more sophisticated random model. However, here we show that our approach with the simplest probabilistic model provides fruitful results.
Below are reported the parts added in the manuscript to answer to this point.
Relevant changes in the Discussion section:
“To design the probabilistic model, we use a minimalistic approach. We use the assumption that directed connections are represented by the matrix of independent Bernoulli random variables. […] One way to overcome this limitation might be the use of more sophisticated probabilistic processes where the random variables corresponding to different connections become dependent (e.g. random Markov field approach).”
4) The last paragraph in the Discussion makes it seem as if the authors are proposing a novel methodology. However, this novelty is not clear since employing probability matrices to define models is a standard practice: its a feature included in several major neural simulators (NEST, Brian, Moose,.…), and many existing models (e.g. on ModelDB) employ prob matrices, either at a populationlevel or celllevel (for small networks and when enough exp data is available).
We agree that the last paragraph was not clearly formulated and have adjusted it. We also agree that specification of connection probabilities is a simple and general idea which was used by several authors and neural simulators for generating neuronal connectivity. The novelty of probabilistic model is that it can be used to compute connection probabilities. In fact, we have three statements here:
1) Tadpole’s spinal cord is a unique case for the study of neural connectivity, which allows us to develop both a biologically realistic connectome and its generalization, the probabilistic model.
2) Generally, finding the biologically realistic matrix of connection probabilities is a very difficult problem. There are several publications attempting to define the probabilities of connection using limited biological data and we provide references to these publications at the last paragraph of Discussion section.
3) We discuss how to use the nonlinear optimization technique to find the connection probabilities under some biological constraints. To adjust probabilities, a cost function can be derived to take into account biological limitations. Optimization of this cost function will provide the best probability values.
We replaced the last paragraph of the Discussion section with the following:
“Finding neural connection probabilities under biological constraints is a difficult problem. […] We will report this result in a separate publication.”
5) Authors should include more details about the functional model employed since this is crucial for understanding the results (e.g. low connection weights can drastically alter the effect of high connection probabilities). A table with the main parameters used in the functional model should be included. Additionally, authors must share all the model and analysis code  this is common practice nowadays in computational neuroscience to ensure replicability and reproducibility.
Details about the anatomical and functional model formulation and parameters were added in the Appendix.
We provided the code for running all the simulations of the anatomical, probabilistic and functional models in ModelDB. We have temporarily made the code private, but we will make it public once the work has been published. Even though the code is private, it can be downloaded using a readonly “access code”. To download the code, use the link and access code given below:
https://senselab.med.yale.edu/ModelDB/enterCode.cshtml?model=238332
Access Code: tadpole1
The folder “spinalcordmaster/functional model/paper figures” contains the files for generating the figures in the paper.
To run these files, it is required to first generate a bunch of anatomical connectomes, the probabilistic matrix and the functional model. To do so:
1) Enter in the tadpolespinalcordmaster/anatomical connectome folder and run the MATLAB file run_many_connectomes.m specifying the number of anatomical connectomes to generate (this will generate output anatomical connectome files is the folder “connectome files”).
2) Run the Python file tadpolespinalcordmaster/probabilistic connectome/build_probabilistic_model.py, specifying the number of anatomical connectomes to use for the probabilistic model. This will automatically generate the probabilistic and anatomical model files in spinalcordmaster/probabilistic connectome/anatomical adjacency matrixes.
3) Enter the folder “spinalcordmaster/functional model/” and compile the nrnmod files in the nrnMod folder (via nrnivmodl Neuron command). The file main.py is the main for launching the simulations and generate output files.
For each of the models a README file contains additional explanations on how to run the simulations
[Editors' note: further revisions were requested prior to acceptance, as described below.]
The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:
1) The explanation of why the probabilistic model is an advance over the anatomical model is not yet convincing. Looking at both the response to the reviewers and the argument in the Introduction, it still seems that the same results could have been obtained by consideration of the distribution of networks generated by the anatomical model. The authors must clarify why they emphasize the need to construct the probabilistic matrix as a new step in the analysis process, otherwise averaging the outcome of an existing model and putting into a simple probabilistic framework is a modest and intuitive extension of the previous work without added value.
We are grateful to all reviewers for these fruitful comments, which clearly show that the motivation and advance of the probabilistic model should be better explained:
The probabilistic model provides a different way to look at the information generated by the anatomical model. It is grounded in the previous anatomical model as the anatomical model is grounded in the biological anatomy. It provides a different perspective on data generated by many anatomical models. Connectivity analysis could be performed directly on results of the anatomical model, however, the anatomical realisations are highly detailed and variable making it actually very difficult to analyse their structure and determine what is important and what is not. In the probabilistic metamodel all connections are independent, again aiding analysis of the generalised structure and allowing identification of core properties. The probabilistic model is based on multiple anatomical networks and simplified by picking up only two elements (the order of neurons along the rostrocaudal dimension and the contact probabilities) to generalise their structural properties. It is the different perspective that makes the probabilistic model an advance.
Also, we believe that the following text highlights the main advances of the probabilistic model and we added this text after the Discussion section as a new “Conclusion” section:
“We study the structure and function of the spinal cord neuronal network using experimental data and computational modelling. Our anatomical model generates multiple highly variable and nonhomogeneous connectomes and to deal with this large and complex data we design a very simple mathematical metamodel expecting that this new probabilistic model will reflect (generalise) structural properties of anatomical connectomes and show proper functioning. […]
et al.,The probabilistic model provides a different way to look at the information generated by the anatomical model. It is grounded in the previous anatomical model as the anatomical model is grounded in the biological anatomy. It provides a different perspective on data generated by many anatomical models, and it is this different perspective that makes the probabilistic model an advance.”
2) The authors have not fully addressed the issue of why the celllevel probabilistic model is capturing general/fundamental structural features of the tadpole network and not overfitting.
See our responses to reviewer #3 comment
Without addressing these two issues more thoroughly, the impact of the paper rides on its scientific results, connecting the structure of the model to key aspects of its functional output and the authors' emphasis on method is misplaced. Some clarifying points from the expert reviews are presented below.
Reviewer #3:
I don't think the authors have addressed the issue of why the celllevel probabilistic model is capturing general/fundamental structural features of the tadpole network and not overfitting. They argue that it is the result of averaging over 1000 anatomical model connectomes, but these are still a model, e.g. they have a fixed number of neurons for each population, which is not the case for real tadpole networks, and so describing the structure at the celllevel prevents it from being applicable to specimens with a different number of neurons. They also argue that all connectomes achieve reliable swimming, but I fail to see how this is relevant. I think the paper still has value despite this limitation, but I think it should be clearly stated in the paper.
1) See our response to reviewing editor comment 1.
2) Fixed number of neurons:
Currently, in the anatomical model, neuron population sizes and their broad rostrocaudal distributions are fixed, based as far as possible on biological estimates, and are symmetrical between the two sides. Everything else varies from one connectome generation to another, again based on biological measurements: the exact rostrocaudal position of each soma, and therefore the start position of each axon; axon outgrowth angle, the exact trajectory followed, and the length of each axon. We have assumed that the main variation in the distribution of connections will result from these differences in the patterns of axon growth rather than differences in soma numbers, which we expect to be relatively similar between individual animals. Our “averaging” of connectomes should therefore have covered the main variability that would be expected between individuals.
In reality, although biological data are limited, total numbers of spinal neurons and the population sizes of individual neuron types do not appear to vary greatly between animals (perhaps ± 1015% at most) at this early stage of development.
Any variation that there is in numbers is small compared to the differences from the previous developmental stage and the following stage, both of which swim in the same way.
To address this comment, we included a short clarification to Appendix 1 (subsection “Branching and commissural axons”) related to the anatomical model.
3) Overfitting:
We confess to being a little unsure about the use of the term “overfitting” here. According to our understanding, overfitting occurs when the number of parameters is large, and the amount of data is not sufficient for parameter estimation, meaning that one set of possible parameters is chosen when in fact many would produce the same fit to the data. In the probabilistic model, the number of parameters is of order 1500*1500 (number of connection probabilities), whereas the amount of data for estimation is about 1500*1500*1000 (number of pairs of neurons across 1,000 anatomical connectomes). The probabilistic model is therefore, according to our understanding, a proper generalization of the data that does not suffer from overfitting. If the reviewer can provide a clearer description of their concern, we will be happy to address it in more detail. Our attempt to address the overfitting matter was clearly not successful and has led to confusion, therefore we deleted the sentence from the Introduction.
4) Probabilistic model is capturing general /fundamental structural features:
We agree that our formulation is confusing. In fact, the fundamental feature is the ability to generate reliable swimming by all probabilistic connectomes. The text (Introduction) has been adjusted.
https://doi.org/10.7554/eLife.33281.019Article and author information
Author details
Funding
Biotechnology and Biological Sciences Research Council (BB/L000814/1, BB/L002353/1)
 Andrea Ferrario
 Stephen R Soffe
 Roman Borisyuk
Plymouth University
 Andrea Ferrario
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Prof Alan Roberts for his valuable contribution to model development, discussion of results and many other useful comments. We also thank Dr WenChang Li for his valuable assistance. This work was supported by the UK Biotechnology and Biological Sciences Research Council (BB/L002353/1; BB/L000814/1). AF was supported by a PhD studentship from Plymouth University. We would like to thank the Reviewing Editor Prof Ronald Calabrese and anonymous reviewers for their fruitful comments. RB is on leave from the Institute of Mathematical Problems of Biology, The Branch of Keldysh Institute of Applied Mathematics of Russian Academy of Sciences.
Reviewing Editor
 Ronald L Calabrese, Emory University, United States
Publication history
 Received: November 2, 2017
 Accepted: March 25, 2018
 Accepted Manuscript published: March 28, 2018 (version 1)
 Version of Record published: April 20, 2018 (version 2)
Copyright
© 2018, Ferrario et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,274
 Page views

 185
 Downloads

 5
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.