Synaptic Strengths Dominate Phasing of Motor Circuit: Intrinsic Conductances of Neuron Types Need Not Vary across Animals

Visual Abstract


Introduction
Through their projections and synapses onto motor neurons, central pattern generator (CPG) circuits control motor output for a variety of rhythmic behaviors. Across individuals of a species, motor circuits, both motor neurons and CPG elements, show large variability, not only in intrinsic neuronal parameters (Bucher et al., 2005), but also in circuit synaptic parameters (Goaillard et al., 2009). Given this animal-to-animal variability, how do motor circuits produce functional, albeit individually unique, output patterns across animals (Norris et al., 2007a;Goaillard et al., 2009;Williams et al., 2013)?
This question has been largely addressed, not with experimental animals coming from controlled backgrounds, but with those collected from their natural habitats (Goaillard et al., 2009), which have multiple factors that can contribute to biological variability. These genetically unique individuals must adapt to developmental and environmental variability presumably by applying some rules (e.g., homeostatic or developmental) to achieve functional output. We are interested in explaining how such rules act at the neuronal and circuit level to achieve their unique functional network output. In particular, what are the relative contributions of synaptic and intrinsic parameters in the production of unique functional output across animals?
Central pattern generating networks of invertebrates have provided some of the best evidence showing functional output with underlying variability of intrinsic and synaptic conductances (Prinz et al., 2004;Goaillard et al., 2009;Ransdell et al., 2013;Marder et al., 2015;Doloc-Mihu and Calabrese, 2016). In CPGs, the firing phase of component neurons is considered a critical aspect of a functional motor pattern, and phase varies considerably, albeit within functional limits, across animals as shown in the stomatogastric (STG) nervous system (STNS) and in the leech heartbeat system (Wenning et al., 2004a,b;Bucher et al., 2005;Norris et al., 2006Norris et al., , 2007b. For leech heartbeat, each animal arrives at a unique solution to produce a functional heartbeat motor pattern based on phase differences in the premotor pattern and synaptic strength patterns from the CPG premotor interneurons to motor neurons (Norris et al., 2011;Wright and Calabrese, 2011b).
Previous large-scale modeling simulations in other systems, which included both synaptic and intrinsic channel variations, generated databases of different conductance parameter configurations, but these were compared to population averages rather than data collected from individuals (Prinz et al., 2004). Databases of simplified computer models were also constructed to explain variable phasing observed in a multianimal dataset (Williams et al., 2013). In this article, we focus on the coordination of heart motor neurons (HEs) that each receive their inputs from inhibitory premotor heart interneurons (HNs). While all input-output spiking phases and synaptic strengths were recorded, intrinsic conductances in these individuals could not be measured simultaneously. We use computer models to match neuronal profiles (i.e., specific biological instances of network spike phasing and synaptic parameters) to reverse engineer motor neuron intrinsic properties, such as in the case for single animals (Lamb and Calabrese, 2013), but extend this methodology to multiple individuals.
Premotor and motor phasing and synaptic inputs to leech heartbeat motor neurons were experimentally characterized from a set of six individuals (Norris et al., 2007a, b;Wright and Calabrese, 2011b). The bilateral motor circuit was reconstructed by using the same parameterized conductance-based multicompartmental motor neuron model as in Lamb and Calabrese (2013). We optimized sets of model instances (with varying parameter values) of the model neuronal circuit using a multiobjective evolutionary algorithm (MOEA) within the observed natural variability to produce each individual leech's unique motor output and assessed the relative contribution of intrinsic ionic conductances versus synaptic parameters. We started by keeping synaptic strengths at the measured average values for the individual and allowing only variation of the intrinsic ionic conductances and bilateral electrical coupling between motor neuron pairs. We found that functional models were difficult to find and that they were not successful when tested with the premotor inputs and expected motor output from other individuals in our set. However, allowing also the synaptic strength parameters for an individual to vary within the error of our measurements for that individual led to much higher success rate for our motor neuron model searches. Moreover, such motor neuron models were highly successful when challenged with premotor parameters (synaptic strengths allowed to vary within the measurements' error) from other individuals. Thus, synaptic strength parameters appear more critical than intrinsic membrane parameters in the expression of animal-to-animal variation in motor output in the leech heartbeat system.

Electrophysiology and spike-triggered averaging (STA) of synaptic currents
Electrophysiological recordings were obtained from previously published data (Norris et al., 2011). Recordings were made in isolated nerve cords of Hirudo sp. extracellularly from premotor HN interneurons and in voltageclamp from ipsilateral segmental HE motor neurons. Norris et al. (2011) described 12 preparations that have complete recordings of input spiking pattern (premotor interneurons), output spiking pattern (motor neurons) and synaptic strengths in two motor neurons of all four premotor interneurons in both peristaltic and synchronous coordinations (see beginning of Results section for their definitions). They showed statistically that the synaptic strengths of these preparations were representative of a sample involving many more preparations of less complete data from previous studies (Norris et al., 2006(Norris et al., , 2007a. Of these 12, six were perfect in the sense that there were no changes in period between the recordings of the two coordinations so that data from one side could be reflected to give a bilateral pattern; this was very important for making a bilateral pattern for the modeling done here. The pattern of interneuron activity and motor neuron activity of these six is representative of the 12 and enabled very detailed modeling and dynamic clamp studies that ultimately led to this study (Wright and Calabrese, 2011a,b;Lamb and Calabrese 2013). A recently published meta-analysis of all the preparations (Ͼ100) is available in a project database with all data available on Dryad (https://datadryad.org/resource/doi:10.5061/dryad. c0g0p/3; Wenning et al., 2018) enabling external validation of the representativeness of this sample. We show here that the premotor patterns and motor patterns in both coordinations are representative of the much larger sample in the project database (Extended Data Fig. 2-1).
Recorded electrophysiological activity traces were analyzed using custom scripts written in the MATLAB computing environment (MathWorks). These analysis scripts extracted bursting characteristics of recorded neurons (https:// github.com/RonCalabreseLab/Burst-Phase-Matlab-lib). They were also applied to quantify characteristics of simulated neurons (see below). These characteristics were used as the target metrics to evaluate whether model neurons are similar to recorded ones. The analysis started by finding beginning and ending spikes in each burst, their median spikes (indicated with diamonds in figures) that is used to compare phasing of bursts, within-burst spike rate and amplitude. It also low-pass filtered the bursting activity to extract the slow-wave envelope on which the bursts rode. The height of this slow-wave envelope was also one of the target metrics.
For assessing the strength and dynamics of synapses, we have employed an offline spike-triggered averaging (STA) method to analyze IPSCs caused by premotor HNs that are recorded in voltage clamp from HEs, similar to a previous method (Norris et al., 2007a). STA was used to determine synaptic connectivity, strength, and dynamics (short term synaptic plasticity) in both coordination modes, but using data recorded from one side of each preparation. Data from the two coordination modes were combined to increase the signal-to-noise ratio in STA because no significant differences were noted between coordination modes in synaptic connectivity, strength, or dynamics (Norris et al., 2007a).
Spike and IPSC detection/averaging were performed off-line using custom MATLAB software (https://github.com/ RonCalabreseLab/HN-HE-synapses-STA). For most of the duration of their bursts, the premotor interneurons tend to fire at a quasi-constant rate, often resulting in multiple peaks in the STA. The largest of these peaks was selected in each trace for measurement of amplitude and latency as before (Norris et al., 2007a). Data used in the STA estimation can be accessed publicly (https://doi.org/ 10.6084/m9.figshare.7963757).
We have improvements over the previously reported STA method in several aspects: (1) we verified that proper number of leading and trailing IPSCs are removed from the STA in each burst; (2) we have an interactive feature to exclude from averaging IPSCs that are near post-synaptic spikes, which escape from voltage clamp because of poor space clamp of the spike initiation zone; (3) we provide analysis of within-animal variability of the STA and its comparison to the average (Fig. 1C); (4) finally, we use an improved method to align the current offset between traces from different bursts. Previous methods reported align traces to have zero current offset at their value that corresponded to the presynaptic spike trigger time of the selected HN. However, the trace could shift up or down considerably with this method based on the coincidental existence of IPSCs from other HNs during this time. Instead of this method, we used a bandpass filter (5-3000 Hz) to remove low-frequency current offsets from traces independent of their value at the spike trigger time. This method reduced the offset visibly and quantitatively allowing more accurate averaging.

MOEA optimization of a conductance-based ensemble model of HE motor neurons
To simulate HE motor neurons, we used a previouslypublished multicompartmental conductance-based model (Lamb and Calabrese, 2013; model files available on the ModelDB repository, https://senselab.med.yale.edu/ModelDB/ ShowModel.cshtml?modelϭ153355). Each motor neuron in this model is composed of seven compartments: a soma, three serially-connected primary neurite compartments, a secondary neurite compartment (passive only) branching from the primary neurite, a synaptic compart-ment connected to the secondary neurite, and an axon serially-connected to the end of the primary neurite. Each of these included combinations of intrinsic ionic channels specified in Table 1. The ensemble model of the heartbeat motor neuron circuit contains bilateral pairs of these motor neuron models coupled together electrically with a conductance (g synE ) at their respective synaptic compartments. Two pairs of model neurons are simulated for leech segments 8 and 12. Different segmental HE motor neurons receive their corresponding recorded HN inputs. We modified this model (https://github.com/RonCalabrese Lab/HE-model) to allow varying passive parameters independently across different compartments, added variable global and individual synaptic multipliers, and allowed user specification of different HN interneuron input patterns recorded from different leeches in our data set (https://doi.org/10.6084/m9.figshare.7963733). We focused our investigations on one model instance from Lamb and Calabrese (2013), whose maximal conductance parameters are shown in Table 2. Lamb and Calabrese (2013) had previously used a MOEA approach to optimizing model circuit parameters to find motor patterns that match recordings from a single animal (5/19B). Model outputs were quantified based on measurements of activity waveforms (Table 3). Target ranges from each animal's target values indicate acceptable metric values during ensemble model fits. We used the same target ranges in this study to maintain consistency and produce comparable results. We have extended this schema to fit model outputs to targets from multiple animals. Target ranges selected can be compared in the table to SDs in peristaltic (peri) and synchronous (sync) coordination mode measurements from HE(8) and HE(12) motor neurons averaged over data recorded from six leeches. Duty cycle and phase metric target ranges were less than one SD over the data set. Frequency had an especially relaxed target range. We ignored the absolute phase metrics and instead used the motor phase progression metric. Motor phase progression depends on the absolute phases of motor neurons because it is the phase difference of HE(8) and HE(12). In total, there were 18 metrics that the optimization algorithm maintained within the target ranges. The analysis was performed using custom MATLAB scripts (https://github.com/RonCalabrese Lab/HE-model-analysis-matlab).
Multiobjective optimization algorithms have been successfully used for conductance model parameter optimization before (Druckmann et al., 2007). The basic idea is to find solutions that optimize each objective separately, instead of previous methods that weigh together different objectives to optimize a summed metric. Optimizing the summed metric may give good solutions in one metric while doing worse in others. The user can specify importance of some metrics by weighing them stronger if known in advance. However, this may skew results in favor of a few metrics. Multiobjective methods behave much more fairly, where solutions are evaluated for their goodness by checking all objectives. Therefore, they are more appropriate for matching biological features, which are all important to match. These methods have become the standard in ambitious neuronal modeling projects (Markram et al., 2015), by accumulating many new advancements (for review, see Druckmann, 2013). One of the most commonly employed MOEA algorithms is the non-dominated sorting genetic algorithm (NSGA-II; Deb et al., 2002). However, this algorithm is known to be unsuitable for problems with more than three objectives (Zou et al., 2008), which makes it a suboptimal choice for our optimization task with 18 output metric objectives. Lamb and Calabrese (2013) have employed a custom modification of NSGA-II that partially eliminated this handicap, the elitist non-dominating vector evaluated genetic algorithm (endVEGA), which was written in the Cϩϩ programming language and used a roulette wheel to choose more of the lower-value objective metrics irrespective of their number (Milanova et al., 2006;Smolinski et al., 2007). Since an easily customizable implementation of this algorithm was not publicly available, we adopted a custom open-source MATLAB toolbox called GODLIKE (http:// www.mathworks.com/matlabcentral/fileexchange/24838godlike-a-robust-single-multi-objective-optimizer). We have customized and improved GODLIKE by allowing parallel processing on multithreaded hardware. GODLIKE offered four multiobjective optimization algorithms that can be combined together, and included an NSGA-II implementation. As expected, it failed to converge in our 18dimensional objective function space, because it could always find solutions where some metric was better, which could not be dominated by any other population member. For example, it could find a model instance that has a perfect match to the target spiking frequency, while several other metrics remained outside of target range. We stopped this overzealous single-objective fitting by preventing the optimization of metrics below a threshold, in particular the limits of our target range for each metric. A fuzzy truncation function on the squared metric value (Extended Data Fig. 3-1A) allowed us to place a lower bound on metric errors equal to the target range value, such that they cannot be improved further. Therefore, the algorithm focuses on improving the remaining metrics (Extended Data Fig. 3-1B). In practice, this method did not dramatically improve the quality of solutions found, but it did improve the speed of convergence and the number of evaluations to reach good quality (low target metric error) solutions (Extended Data Fig. 3-1C,D). That is, many more low maximal absolute error (MAE) solutions were found with same number of evaluations. Control of simulations, comparison and sorting of metrics, and plotting were achieved using the Pandora toolbox (Günay et al., 2009; RRID: SCR_001831) for MATLAB. Figures in this article can be generated using the public resources for the source code (https://github.com/RonCalabreseLab/HEmodel-animals-figures) and the associated simulation data (https://doi.org/10.6084/m9.figshare.7963748).

Statistics
Synaptic strengths measured from living animals and computer simulations were compared using the Kruskal-Wallis test. Being a non-parametric test, Kruskal-Wallis does not assume samples follow a specific distribution. It tests the null hypothesis that medians between the two distributions are equivalent. We chose this test as it is more stringent than a one-way ANOVA test, which assumes the underlying distributions are normal. The distributions of the living and simulated synaptic values were shown in detail in a supplementary figure (Ex-tended Data Fig. 6-2). The Kruskal-Wallis test was applied using the implementation found in MATLAB (MathWorks) named "kruskalwallis." Correlations between intrinsic P and K 2 conductances were tested using multiple linear regression in MATLAB (MathWorks) with the "regress" command. R 2 A B C D Figure 1. Experimental data from multiple animals show large variability in interneuron phase, output phase, and synaptic strengths. A, Representative ipsilateral connectivity schematic and extracellular traces from preparation 5/27B in peristaltic coordination showing premotor heart interneuron (HN) neurons in segments 3, 4, 6, and 7, with projections to heart motor neurons (HEs) in segments 8 and 12. Phase difference in the premotor inhibitory synaptic input (pm⌬) from heart interneuron HN(x) to HE(y) results in a motor phase progression (m⌬ x,y ). B, Phases of both HNs and HEs vary across preparations. C, STAs of relative synaptic strengths of HN(3) and HN (7) to HE(8) synapses in preparation 5/27B showing the mean (black lines) and SD (gray areas). Red stars show the points between which amplitudes are measured. These two synaptic currents correspond to the red stars in panel D. D, HN to HE relative synaptic strengths (mean Ϯ SD) from six individual animals show variability in their synaptic strengths and disparity between inputs and their ordering of magnitudes. and p statistics were reported from the output of this command.

Experimental data from multiple animals vary in interneuron phase, output phase, and synaptic strengths
The leech heartbeat is a bilateral system, where blood is pumped by rear-to-front peristalsis on one side (heart tube) and nearly synchronously on the other side (heart tube). This activity pattern of two distinct coordinations, peristaltic and synchronous, switches sides periodically every few minutes to ensure that blood circulates equally on the two sides. The same bilaterally paired premotor interneurons and motor neurons produce both coordinations, and in particular only the relative timing of interneurons and motor neurons changes between the two coordinations ( Fig. 1A; Norris et al., 2007aNorris et al., ,b, 2011. It has recently been shown that there are significant side-to-side differences in phasing of activity for each coordination and in underlying synaptic strengths but on a given side, motor output is very similar across repetitions of a coordination (Wenning et al., 2018).
We focus on the coordination of motor neurons (HEs) in segments 8 and 12 that each receive all of their inputs from the same four inhibitory premotor interneurons (HNs) in segments 3, 4, 6, and 7. When one side is in peristaltic coordination, bursting spike patterns of premotor interneu-rons [HN(3)-HN (7)] show a rear-to-front phase progression (Fig. 1A). Bursting phase of premotor interneurons is highly variable across animals and controls the phasing of motor neuron [HE(8) and HE (12)] output spiking patterns (Norris et al., 2007a). In some animals (e.g., in 5/19B), even premotor phase ordering may be different (Fig. 1B), but the general trend of rear to front progression is always observed. Variability in premotor input phasing is also reflected by varying motor phases of HE(8) and HE(12) in each animal, but the premotor phase difference (pm⌬) always sets an upper limit on the forward motor phase progression (m⌬; see Wright and Calabrese, 2011b). This creates the rear-to-front peristaltic wave constricting heart muscle that is critical for pumping blood and therefore the survival of the leech. Functionality of this feature, despite differing input patterns, depends on the strengths of synapses connecting the interneurons to motor neurons, which also varied across animals.
Previously, Norris et al. (2007a) estimated and reported synaptic strengths using a STA method. Here, we are also reporting within-animal variability of synaptic strengths ( Fig. 1C) using an improved STA method (see Materials and Methods). Synaptic strengths were normalized per animal to have a meaningful comparison across animals and preparations. When looking across animals, relative strengths were, not only highly variable, but also ordered inconsistently by their magnitudes (Fig.  1D). Motor neuron HE(8) had a more even distribution of relative synaptic strengths across all four inputs compared to HE(12) that had small versus large strengths. Across animals, the strongest input to HE(8) alternated between HN(4) and HN(6). Inputs to HE(12) were more consistent across animals, HN(6) and HN(7) being always the strongest; however, their ratios still varied greatly. Other than these rules, synaptic strength patterns seemed to be unique for each animal despite the criticality of their function.

Functional metrics of the heartbeat circuit such as the motor phase progression are also variable
Functional output metrics of the heartbeat circuit also showed variability across animals in both coordinations ( Fig. 2; for a comparison of these distributions to the larger project dataset, see Extended Data Fig. 2-1). These metrics were previously described in detail (Lamb and   These ranges provided the limits around each animal's target values for fitting ensemble models. Phase and duty cycle are displayed in normalized units. Metrics of slow wave height and spike height were used to ensure models were physiologic and were estimated from recordings. Target ranges are compared to mean and SD of the measurements across six animals for HE(8) and HE(12) in the columns to the right for peristaltic (peri) and synchronous (sync) coordination modes. Motor phase progression is calculated using phase of both motor neurons, so it was displayed only once under HE(8).
Calabrese, 2013), and include duty cycle, firing rate, slow wave amplitude, and spike rate. The most important output metric of a functional heartbeat is the peristaltic motor phase progression (m⌬), which was calculated here as the phase difference from motor neuron HE(8) to HE(12) and exhibited substantial animal-to-animal variability. Variability also existed in the synchronous coordination, but the phase differences were much smaller. Bursting duty cycle and firing rate of the motor neurons were similarly variable across animals. Each animal having different input phase progression, synaptic strengths, and output phase progression, raises the question whether output variability is caused by inputs only, or intrinsic ionic conductances must also differ from animal to animal. Another CPG circuit, the STNS also exhibits similarly stereotypical but variable neuronal output activity in the same neuronal types across individual crustaceans (Bucher et al., 2005), while underlying channel conductances and gene expression patterns vary significantly Schulz et al., 2006). The two to four-fold difference found in maximal conductance magnitudes has been suggested to cause the observed population variability in neuronal activity. Variable underlying conductances in neuronal types was also observed in mammalian Purkinje cells (Swensen and Bean, 2005). We can investigate how changes in inputs and ionic conductances contribute to explaining the outputs in the leech by using a computer model of the heartbeat motor circuit.
No model instance produced functional output using the measured synaptic weight averages Previously, Lamb and Calabrese (2013) constructed a conductance-based multicompartmental model of the heart motor neurons (HEs). This is an ensemble model that simulates both sides of the animal because bilateral HE motor neurons are electrically coupled. Premotor input phase and synaptic strength recorded from one side (for both activity coordinations) in experiments were used as inputs to both sides of the model circuit. Synaptic strengths were previously measured by STA experimentally for each HE motor neuron independently (Wright and Calabrese, 2011b). Not knowing how strengths of different HE synapses compare to each other across animals, experimental strengths were interpreted in relative terms. Here, relative synaptic strengths to each HE were scaled by a global scaling parameter () to generate absolute strengths to be used in simulations (Fig. 3A). The model could then be simulated by playing back inputs from model synaptic currents triggered by recorded HN spike patterns, aiming to recreate HE output phasing that matched recorded animal's synaptic profiles. This simulation method allows us to explain how the measured properties of the circuit come together to generate motor output and whether it can reproduce observed functional metrics such as burst phase.
Using this method, Lamb and Calabrese (2013) were able to simulate a small motor phase progression target preparation 5/19B (Fig. 1B, left-facing triangles) by using measured inputs to the model (Fig. 3B). However, the model failed to reproduce larger phase progression targets, such as from 5/22B (Fig. 1B, diamonds), using its corresponding inputs (Fig. 3B; compare phase progression values in Fig. 2). Theoretically, the model could achieve large motor phase progression when we artificially maximized output phase progression by strengthening earliest (in phase) and latest synaptic inputs to HE(12) and HE (8), respectively (Extended Data Fig. 3-2A).
Using this method, we confirmed that the output motor phase progression is always limited by, and follows, the input premotor phase progression (Extended Data Fig.  3-2B). However, while the model could match larger motor phase progression values, it was unable to keep other functional metrics within target range for any combination of values for the two synaptic scaling parameters () of HE(8) and HE (12), respectively (Extended Data Fig. 3-2C). Furthermore, searching across all maximal conductances (for ionic currents and bilateral electrical coupling between HEs) and HE synaptic scaling parameters of the model failed to find any solutions within all metric target ranges for preparation 5/22B (Extended Data Fig. 3-3). Since intrinsic ionic conductances could not compensate for the model's shortcomings, we turned back to modifying recorded relative synaptic weights from their average values measured by STA.

Slightly adjusted synaptic weights result in functional models without changing intrinsic ionic conductances
By preserving the same intrinsic ionic maximal conductance magnitudes in the model instance and allowing synaptic strengths to slightly vary from their measured averages, we were able to find functional models for target 5/22B (Fig. 4) and 5/20B (Fig. 5). To achieve this, we defined synaptic multiplier parameters from each HN input to an HE motor neuron (Fig. 4A) and allowed them to vary only within Ϯ1 SD of the recorded synaptic strength averages (Fig. 1D). The synaptic strengths found in the top functional models were within Ϯ1 SD of the measured averages for preparation 5/22B (Fig. 4B). The simulated models were able to produce the target phase progression of the individual animal with no other electrophysiological anomalies (Fig. 4C). The measured functional metrics were all within target ranges as well (Fig. 4D). If we had used the original measured synaptic strength averages, the metrics would be outside of this range for the synchronous duty cycle of HE(12) and much lower than the target of the peristaltic motor phase progression (Fig.  4E). To show another example, we were able to find synaptic strengths (Fig. 5A) and functional output metrics (Fig. 5C) within Ϯ1 SD of the measured targets for preparation 5/20B that produced proper bursting patterns (Fig.  5B). These two examples gave us confidence for searching nearby synaptic strengths to those measured in the rest of the animals.

Adjusting synaptic strengths resulted in functional models for all six preparations
We repeated the procedure for the remaining four preparations and obtained functional models by finding new synaptic strengths within Ϯ1 SD of measured averages (Fig. 6). The best motor phase progression errors obtained were Ͻ0.6 of its target range across all preparations and largest error was associated with different metrics in each one (Table 4). Bursting activity was also regular across all preparations (Extended Data Fig. 6-1). When compared across six animals, we found no significant difference between relative synaptic strength distributions of measured averages versus topmost models (Kruskal-Wallis, p Ͼ 0.05; Fig. 6A). Furthermore, the model synaptic weights followed the same relative pattern (same distribution) across eight inputs as in living HE neurons (for synaptic strengths found for individual animals, see Extended Data Fig. 6-2). To explain the significance of the synaptic weight pattern, we used the synaptic strength index (SSI) defined by Wright and Calabrese (2011b). The SSI metric calculates the sum of preferences of HN(4) to HE(8) and of HN(7) to HE(12) in a single animal by using where w are relative synaptic strengths. When we applied this equation to both measured and model synaptic strengths for all preparations, we found that model SSI values increased for preparations with larger phase progression, while measured SSI values remained almost flat ( Fig. 6B; see Discussion). Despite this discrepancy and having reproduced the observed peristaltic phase progression values in the model without sacrificing any other functional metric, we turned to the question of whether intrinsic ionic conductances also contributes to individual variability across animals. Surprisingly, all the results we presented so far (except Extended Data Fig. 3-3) were obtained with the same combination of ionic conduc-continued Data Fig. 3-1). Sample model instance (#47) shows bursting patterns where diamonds mark the middle spike of the bursts and the gray triangle connecting the traces of HE(8) to HE(12) indicate the motor phase progression. Peristaltic and synchronous coordination modes are indicated with magenta and blue traces throughout this article. The horizontal sequence of dots atop traces aligned with HN interneuron labels indicate their firing times, while the dots' thicknesses indicate the relative strength between the HN input onto the corresponding HE motor neuron. C, When -scaled ( HE8 ϭ 0.73 and HE12 ϭ 0.75) HN firing patterns of a different animal (5/22B, with target m⌬ ϭ 0.105 Ϯ 0.06) were simulated as input in the HE circuit model, it resulted in a smaller m⌬ ϭ 0.008, which failed to reproduce the increased target value. Extended Data Figure 3-2 gives more details about the impact of the factors on the activity characteristics. Furthermore, searching across all maximal conductances (of intrinsic ionic currents and bilateral electrical coupling between HEs) and HE synaptic scaling parameters of the model failed to find any solutions within one target range for all metrics for preparation 5/22B (Extended Data Fig. 3-3). tances in the model. This suggests that synaptic strengths dominate the phasing of the output and ionic conductances need not vary across animals. Next, we investigated which combinations of ionic conductances resulted in functional models while keeping the newly found synaptic strengths constant.

Intrinsic conductances vary smoothly around synaptic solutions found
So far, we have tested all animal targets with models that employed the same one configuration of intrinsic conductances. To test the impact of intrinsic conductance deviation from this configuration in the heart motor neuron output across animals, we varied intrinsic conductances without allowing changes in synaptic strengths from the newly found values. As expected, we found ensembles of model solutions for each of the animals (Table 5). For each animal, we evaluated at least 5000 model instances in randomized MOEA searches. Model instances within target range ("good" models) were just a small percentage of the total number of instances evaluated, varying between 1% and 13%. This reinforced the understanding that intrinsic conductances are important for producing functional network activity. However, we wanted to assess the scale of this importance and whether intrinsic conductances must vary across animals. Therefore, we asked whether there are other fixed sets of intrinsic conductances that can produce model instances that function across animals. We selected animal 5/22B, which was the most difficult phase target to achieve and which had the lowest ratio of good models. We tested each of its 122 good model instances against other animal targets by setting synaptic strengths to the previously found values for that animal and keeping them constant. . Target 5/22B solution was found by manipulating individual synaptic weights, while keeping the same intrinsic maximal conductance magnitudes of the model instance. A, Individual synaptic multipliers () allow adjusting strengths around recorded averages. B, For input dataset 5/22B, we found 2088 model instances with synaptic weights that are Ͻ1 within-animal SD of target weights while keeping all other functional metrics also within recorded target ranges. C, Bursting activity of models did not show any abnormalities (example model at rank #4 is shown). D, Metrics of same model compared to target. Error bars show Ϯ1 SD across all the targets. Since synchronous phase progression was small and always within target range, we only show peristaltic progression. E, Same model instance with unmodified target synaptic weights was unable to achieve target values. Figure shows best possible metrics achieved by the model when the global synaptic multipliers were scaled to HE8 ϭ 0.7961 and HE12 ϭ 0.7674.
At least some of these instances produced network activity patterns within target ranges for all animals (Table 6, best MAE varied between 0.51 and 0.70). And Ͼ50% of these instances worked in all other animals except target 5/27B (only 28% of 5/22B good model instances were also good for target 5/27B). Having found several intrinsic conductance sets that work across animals supported our view of their secondary role in individual variation of network activity.
Next, we looked at correlation properties across intrinsic conductances in the good sets of model instances that we had found for each animal independently (Table 5). Hudson and Prinz (2010) had suggested that ionic conductance ratios may represent neuronal type and must remain constant across individuals of a population. Here, we provide support for this hypothesis by finding a similar positive correlation pattern across six preparations (Fig. 7) for one of the known ionic conductance pairs (delayed-A B C Figure 5. Adjusted synapses for target 5/20B for the model instance that has the same intrinsic maximal conductance magnitudes as in Figure 4. A, Relative synapse strengths (mean Ϯ SD) recorded from preparation 5/20B (green bars) compared to 3717 models (mean Ϯ SD, white bars) found after performing 10,000 MOEA model evaluations. Traces (B) and relevant metrics (C) of the third highest ranking model. Figure 6. Summary of model synaptic strengths found for all targets. Bursting activity was regular across all preparations (Extended Data Fig. 6-1). A, Boxplots compare distributions of synaptic strengths of recorded (Exps.) versus best performing model (Models) across all six animals. Extended Data Figure 6-2 shows synaptic strengths found for individual animals. B, Synaptic strength index (SSI) compared across targets and mean Ϯ SE of top 10 models (black boxes).

A B
rectifier K 2 and persistent sodium P current in the neurite compartments) that was previously shown to be correlated (Lamb and Calabrese, 2013). We looked for intrinsic conductance correlations in the same data that we presented above (Tables 5, 6). We have found other intrinsic conductance correlations in our dataset, but focused on one conductance pair to perform comparisons across animals.

Discussion
One of the principal questions in neuroscience is how reliable functional neuronal output can be produced while underlying neuronal and synaptic properties are variable across members of a species (Hamood and Marder, 2014). CPG neuronal networks have provided excellent opportunities to study this question (Marder and Calabrese, 1996). Studies across animals focused on properties such as constancy of phase (Hooper, 1997;Bucher et al., 2005;Goaillard et al., 2009;Hamood et al., 2015), period regulation (Kushinsky et al., 2019), and morphologic variability (Otopalik et al., 2017) in the pyloric and cardiac motor patterns of the STG of crustaceans. Similar studies have focused on swimming motor pattern in the lamprey (Grillner, 1974) and on the crayfish swimmeret system (Smarandache et al., 2009). In the leech heartbeat system, the functional output is the phase difference between motor neuron discharge (bursts) in different segments, which control blood circulation through heart tubes. Within functional limits, the phase of the activity patterns of the premotor interneurons and motor neurons, and the synaptic strengths between them, are variable across individuals (Norris et al., 2011), the effects of which are not well understood. Computer models give us the ability to combine disparate experimental measurements to simulate circuit function, which can be used to explore these interactions (Wright and Calabrese, 2011b). We set the parameters of a heart motor neuron circuit model, controlled by the heartbeat CPG, to data collected from a set of six individuals. Our results showed that synaptic strengths play a critical role in producing highly variable animal-specific motor output burst phasing. This is not surprising as the inhibitory synapses from premotor CPG interneurons are responsible for shaping the phasing of the output (Wright and Calabrese, 2011a,b). This is consistent with other studies that look at the interplay between synapses and intrinsic conductances in living neurons across animals (Williams et al., 2013;Anwar et al., 2015;Lane et al., 2016) and across instances of neuronal circuit models (Prinz et al., 2004;Doloc-Mihu and Calabrese, 2011;Williams et al., 2013;Lane et al., 2016). In particular, Lane et al. (2016) showed that, in the crustacean cardiac ganglion, removing one intrinsic conductance can be compensated by changes in another intrinsic conductance and in strength of electrical coupling. When we tested the opposite relation in our model, modification of synaptic strengths could not be compensated by changes in intrinsic membrane conductances.
In our neuronal network model, the values of synaptic strengths measured (as averages in STA) previously did not produce functional output, but we were able to find nearby values within the error of the STA that did. Once functional synaptic weights were found, ionic conductances could be varied to generate ensemble model sets.
Outside of these solution sets, models produce nonfunctional output even with the functional synaptic strengths. Thus, while synaptic strengths predominate in  (12) Best m⌬ value and its error are from the best out of top five lowest MAE (MAE in any metric) models. All top models were within target range error for all metrics. MAE column is for the top model. To avoid overfitting, searches were divided into separate batches, each starting from random initial conditions. For each animal, five batches of searches were executed with 1000 model instance evaluations in each. r g is the ratio of good (functional) models to all evaluated models. ‫22/5ء‬B was searched with five batches of 2000 evaluations. Linear regression was applied to check for correlation between K 2 and P channel conductances, for which statistics are presented (R 2 and p values). Range of metric MAE shows closeness of activity to target range (Ͻ1 accepted as good), "# Good" shows number of functional models, and ratio is the percentage of good models out of 122.
determining output phasing, membrane ionic conductances play an important role in tuning and translating that synaptic input into functional coordinated neuronal output. These results potentially apply to many other motor networks where CPGs exert strong synaptic control over motor neurons.
Why did using measured synaptic strengths fail to generate functional circuit models? Since both activity patterns and synaptic strengths were recorded from the same individuals, it was surprising to find that network model outputs did not match measured features when we used the measured STAs of the strengths. There are several possible explanations for this failure.
Using average values from disparate measurements is known to cause some errors if values come from concave spatial distributions (Goldman et al., 2001;Golowasch et al., 2002). The STA is a well-established method that produces a single averaged value (in this case, synaptic strength) and therefore is not subject to this problem (Ito, 2015). However, if individual EPSCs averaged in the STA come from a bimodal distribution, the resulting averages may be unphysiological. We had already improved the accuracy and increased the signal-to-noise ratio of the STA method employed here (see Materials and Methods), although there could be remaining sources of variability in terms of synaptic plasticity or baseline shifts that may be causing pernicious errors in the measurements. A small error in parameter values can unbalance a whole model particularly in a situation like that here where it is the relative strengths of inputs that is key; a misestimate in one input can throw off the relative strengths.
Another explanation could be a missing feature in the model, although this is unlikely based on previous scru-tiny. The heart motor neuron circuit model is constructed with a high level of biological detail based on physiologic properties measured specifically from the leech. It contains ten ionic currents characterized from voltage clamp recordings; multiple electrical compartments to represent important anatomic distinctions and spatial separation for the axon, soma, and neurites; chemical synapses that feature short-term plasticity, electrical coupling between the heart motor neurons, and conduction delays for intersegmental spikes arriving from premotor interneurons (García et al., 2008;Wright and Calabrese, 2011b;Lamb and Calabrese, 2013).

The same configuration of model intrinsic conductances can create output patterns of multiple animals when synaptic strengths are near measured (STA) values
The same set of intrinsic ionic conductances used in an instance of the HE motor neuron ensemble model was able to reproduce the individually varying output phasing measured from each of the six individual leeches. Usually, activity in identified neuron types are stereotypical across individuals, but with individual variation. This is in contrast to observed 2-to 4-fold variability of intrinsic ionic conductances found in neurons of one type across animals (Bucher et al., 2005;Goaillard et al., 2009). In both of these studies, animals were collected from a natural population that presumably contained significant biological variability. Thus, it is not surprising to observe variability also in neuronal properties, and it can be attributed to environmental, genetic, or developmental factors. However, from a functional network output perspective, is it necessary for these individuals to have different intrinsic ionic conductance magnitudes? Could one transplant the same neuron into another individual and expect to match . Linear regression fits (black lines) were similar across all preparations for functional models (burgundy dots) among a background of 5000 evaluated maximal intrinsic ionic conductance combinations (yellow dots) between neurite K 2 and P channel conductances. Regression resulted in R 2 Ͼ 0.5 for all preparations, shown in each panel's title. All p values were Ͻ0.01 (Table 5). their individual functional output? Since we cannot answer this question easily experimentally, the closest we can get is through computer simulation.
In addition to showing that one set of intrinsic conductances can produce functional network activity across animals, we also showed that one can find multiple such intrinsic conductance configurations. To achieve this, we evaluated thousands of models in the intrinsic conductance space. Only a small percentage of these models produced functional models (Table 5). This confirmed that intrinsic conductances are also important in generating functional network output. We further showed that some of these intrinsic conductance configurations can produce the observed output in multiple animals. Not all of the configurations we selected from one animal produced the observed (target) activity in all animals, possibly because they were selected randomly and not by specific criteria. These results strongly suggest that one can find a large number of model intrinsic conductance configurations that can produce the individually unique functional activity of multiple animals.
Achieving functional output despite observed variability of neuronal intrinsic ionic conductances in wild-caught animals, for example in neurons of the crab STNS, have been attributed to co-regulation rules for maximal conductances based on electrophysiological recordings (Golowasch et al., 1999a,b) and mRNA counts (Schulz et al., 2006(Schulz et al., , 2007. Co-regulated intrinsic conductances found in neuron types thus often compensate for each other and produce stereotypical activity patterns (Hamood et al., 2015;Lane et al., 2016). Therefore, although intrinsic conductances are variable, output produced is similar, which supports our finding that a single type of neural activity produced by the same set of intrinsic model parameters could function in multiple animals as long as the relative synaptic strengths are appropriate to the individual. Within error limits, the relative synaptic strengths could also vary in a large number of model solutions. By definition, relative synaptic strengths imply that, across neuronal models, absolute values are linearly related. Therefore, the synaptic strengths that we found using the model are also correlated. Having already confirmed experimentally some of these predictions (Norris et al., 2011), our model results can be taken to further analyses.
In a computational modeling study of the crustacean STG, Hudson and Prinz (2010) hypothesized that characteristic neuronal types must have similar co-regulation rules. However, this conclusion was reached by constructing ensemble neuronal models that match population targets rather than data from individual animals. Ensemble neuronal modeling creates a large set of neuronal models that fit a certain output criterion and mine it for biological questions (Prinz, 2010;Doloc-Mihu and Calabrese, 2011;Günay, 2015). Here, we provide support for the hypothesis of Hudson and Prinz (2010) by constructing ensemble model databases of the heartbeat motor neuron circuit model across animals. Once we established the new individually appropriate synaptic strengths, we kept them constant and varied intrinsic ionic conductances to find sets able to match each animal target. In the resulting ensemble datasets, we showed that some intrinsic ionic conductances exhibit similar coregulation rules across animals, supporting the earlier prediction and suggesting a required motor neuron type for functional motor output.
In summary, we find that synaptic strengths are critical in creating the proper output phasing in a rhythmic motor circuit. Rhythmic circuits are responsible for many types of motor activity such as walking, chewing, breeding, and breathing across many species and are critical for life (Marder and Calabrese, 1996). The type of CPG motor circuit that we studied here is similar to circuits in other species and our results may extend beyond the leech. Understanding parameters and limits of animal-to-animal variability would help explain and compensate for artifacts when consistent measurements are required in higherorder animals such as rats (Falcone et al., 2018). Our conclusion is that animals in the wild seem to cope with neuronal and network variability by achieving custom solutions to reach functional output. This is consistent with other studies that show that observed individual variability may not have a significant functional impact (Otopalik et al., 2017). As a result, each individual may have unique combination of synaptic parameters that enable signature stereotypical behavior they express.