Homeostatic Activity-Dependent Tuning of Recurrent Networks for Robust Propagation of Activity

Developing neuronal networks display spontaneous bursts of action potentials that are necessary for circuit organization and tuning. While spontaneous activity has been shown to instruct map formation in sensory circuits, it is unknown whether it plays a role in the organization of motor networks that produce rhythmic output. Using computational modeling, we investigate how recurrent networks of excitatory and inhibitory neuronal populations assemble to produce robust patterns of unidirectional and precisely timed propagating activity during organism locomotion. One example is provided by the motor network in Drosophila larvae, which generates propagating peristaltic waves of muscle contractions during crawling. We examine two activity-dependent models, which tune weak network connectivity based on spontaneous activity patterns: a Hebbian model, where coincident activity in neighboring populations strengthens connections between them; and a homeostatic model, where connections are homeostatically regulated to maintain a constant level of excitatory activity based on spontaneous input. The homeostatic model successfully tunes network connectivity to generate robust activity patterns with appropriate timing relationships between neighboring populations. These timing relationships can be modulated by the properties of spontaneous activity, suggesting its instructive role for generating functional variability in network output. In contrast, the Hebbian model fails to produce the tight timing relationships between neighboring populations required for unidirectional activity propagation, even when additional assumptions are imposed to constrain synaptic growth. These results argue that homeostatic mechanisms are more likely than Hebbian mechanisms to tune weak connectivity based on spontaneous input in a recurrent network for rhythm generation and robust activity propagation. SIGNIFICANCE STATEMENT How are neural circuits organized and tuned to maintain stable function and produce robust output? This task is especially difficult during development, when circuit properties change in response to variable environments and internal states. Many developing circuits exhibit spontaneous activity, but its role in the synaptic organization of motor networks that produce rhythmic output is unknown. We studied a model motor network, that when appropriately tuned, generates propagating activity as during crawling in Drosophila larvae. Based on experimental evidence of activity-dependent tuning of connectivity, we examined plausible mechanisms by which appropriate connectivity emerges. Our results suggest that activity-dependent homeostatic mechanisms are better suited than Hebbian mechanisms for organizing motor network connectivity, and highlight an important difference from sensory areas.


Introduction
Activity-dependent adjustment of nascent synaptic connectivity is a widespread mechanism for tuning network properties in various parts of the developing brain (Kirkby et al., 2013). This activity is relayed by ordered topographic maps to higher processing centers, where it triggers correlated synaptic release from neighboring axonal terminals (Katz and Shatz, 1996;Tritsch et al., 2007;Huberman et al., 2008;Feller, 2009). This timing information is evaluated by postsynaptic neurons to sharpen their connectivity profile through Hebbian plasticity and thus improves network function during this early phase of circuit development.
Locomotor circuits also generate spontaneous activity during embryogenesis in many organisms, including chick, mouse, fish, and fly Hanson and Landmesser, 2004;Crisp et al., 2008;Kirkby et al., 2013). These circuits are not topographically organized, yet endogenous activity patterns and synaptic release are required for coordinated motor behavior to develop (Crisp et al., 2011;Warp et al., 2012;Giachello and Baines, 2015). In Drosophila, motor neurons receive variable amounts of presynaptic input, suggesting that connectivity is plastic and acquires functional patterns in a homeostatically regulated manner (Tripodi et al., 2008;Couton et al., 2015). The underlying mechanisms that might evaluate this early spontaneous activity and establish appropriate neuronal connectivity are unclear.
Computational modeling offers an efficient approach to explore the nature of the mechanisms that could underlie activitydependent tuning of connectivity in a generic motor network that achieves a particular function. Here we study the organization of one such motor network that generates rhythmic behavior, the ventral cord locomotor network of Drosophila larvae. Larvae crawl by peristaltic waves of muscle contractions, which propagate along their body axis (Berni et al., 2012;Heckscher et al., 2012). These coordinated waves travel from posterior to anterior during forward locomotion, and from anterior to posterior during backward locomotion. The locomotor behavior is produced by central pattern generators, which are segmentally repeated and modulated by sensory feedback, but can produce coordinated output even when sensory neurons are removed (Suster and Bate, 2002;Fox et al., 2006;Crisp et al., 2008). Mature crawling behavior of Drosophila larvae is preceded by uncoordinated, yet neurally controlled, muscle contractions and incomplete peristaltic waves during late embryogenesis, which gradually improve coordination before hatching (Suster and Bate, 2002;Crisp et al., 2008). Manipulating endogenous activity during this period can significantly affect the output of the circuit, suggesting that this activity helps to refine connectivity (Crisp et al., 2011;Giachello and Baines, 2015).
We previously established a model for wave propagation, which produces propagating activity patterns with appropriate timing as during crawling in Drosophila larvae (Gjorgjieva et al., 2013). Although distinct connectivity parameters were shown to produce propagating waves of activity, the solution regions are complex (Gjorgjieva et al., 2013;their Fig. 5), and it is unclear how they arise during development. Here, we examine two different activity-dependent mechanisms, which can tune network connectivity during development: (1) a Hebbian plasticity model, where coincident activity between neighboring neuronal populations strengthens connections between them; and (2) a homeostatic model, where synaptic connections are modified to maintain a constant level of excitatory postsynaptic activity based on spontaneous input. We show that homeostatic mechanisms are more appropriate than Hebbian mechanisms to organize connectivity in these motor networks. We demonstrate the case by comparing peristaltic wave properties for the two models based on the final configuration of network connectivity produced by the models. We also demonstrate the robustness of activity-dependent tuning of connectivity by varying the properties of spontaneous activity to predict how manipulations of this activity during embryogenesis might impact output of the ma-ture network. Thus, our work highlights the relative importance of homeostatic mechanisms for establishing functional connectivity in developing motor networks, in contrast to sensory networks where connectivity is more tightly constrained enabling the formation of accurate sensory maps.

Materials and Methods
Network model. The network model was based on the interaction of recurrently coupled excitatory and inhibitory populations to produce robust propagation of activity with appropriate timing relationships as observed experimentally during crawling in Drosophila larvae (Gjorgjieva et al., 2013). Such unidirectional propagation of activity is also observed in other experimental systems (see Discussion). Our model is an abstraction of neural circuits that might generate larval crawling and captures the observed coordination at the segmental level. The model contains eight coupled units of excitatory and inhibitory populations to represent the activity in the eight abdominal segments of the larva (see Fig. 1A).
The activity in each segment was modeled with a Wilson-Cowan unit (Wilson and Cowan, 1972), consisting of two neuronal populations: excitatory ( E) and inhibitory ( I). These two populations represent the joint activity of all central neurons in the central pattern generator circuit for crawling but are sufficiently general to represent excitatory and inhibitory neurons in other systems. The differential equations for the timedependent variation of averaged excitatory and inhibitory neuronal activities are as follows: where the functions G E and G I represent sigmoidal response functions of the excitatory and inhibitory neuronal populations given by the following: where represents the maximum slope of the sigmoid (or if G represents an activation function, it denotes the speed of activation) and represents the location of the maximum slope (or the threshold for activation) ( Table 1). This function was chosen so that G(0) ϭ 0 ensuring that E ϭ 0 and I ϭ 0 is a steady-state solution in the absence of external input (Wilson and Cowan, 1972). The terms k E and k I denote the maxima of the response functions for the excitatory and inhibitory populations ( Table  1). The time constants E and I control decay of excitatory and inhibitory activities after stimulation, and determine the timescale of network activity. The connectivity coefficients a (or e) and c (or f ) represent the average number of excitatory and inhibitory synapses per cell in the excitatory (or inhibitory) population, respectively. The time-varying function P ext (t) denotes the external input applied to an excitatory population. To test whether the mature networks have sufficiently strong connections between segments to propagate a wave, waves were initiated by providing constant input P ext ϭ 1.7 to the most posterior (anterior) segment for forward (backward) waves. During development in both Hebbian and homeostatic models, P ext (t) was a time-dependent function representing the spontaneous input into each segment (see below).
To model wave propagation, we coupled eight Wilson-Cowan E-I units. Two types of connections were created between neighboring segments. First, bidirectional excitatory connections (b) between excitatory populations of neighboring segments allow activity to propagate along the segments. Second, inhibitory connections (d) from the inhibitory population in one segment to the neighboring (anterior and posterior) excitatory populations terminate activity in each previously active segment and ensure unidirectional wave propagation (see Fig. 1). This network was equipped with symmetric bidirectional connectivity such that forward and backward propagating waves were generated with similar properties, consistent with experiments (Gjorgjieva et al., 2013). We list example parameter values that generate a functional wave in Table 1. A "contraction" in the model was defined as suprathreshold activity (activity above a threshold, 0.2 or 0.3) of the excitatory population (see Fig. 1B, inset, horizontal dashed line). Such contractions were used to implement the activity-dependent weight refinement in the Hebbian and homeostatic development models. We examined the timing relationships between neighboring segments during propagating waves of suprathreshold activity in the model by analyzing two quantities: The duty cycle was the total time that each excitatory activity is suprathreshold relative to the wave duration. The interburst interval was the time from the initiation of suprathreshold activity of one segment to the initiation of suprathreshold activity of the neighboring segment, normalized by wave duration.
Weights are labeled with a superscript "f" or "b" to denote their respective role in forward or backward wave propagation, but they are both present in the network at all times (see Fig. 1). Thus, there were 28 modifiable weights: The model was initialized with weak connectivity incapable of generating waves. We considered two scenarios: (1) All the weights started identically b ϭ Ϫd ϭ 2; this was 10% of the values used to generate robust waves (Gjorgjieva et al., 2013). (2) The weights were independent random samples drawn from a uniform distribution in the range between 0 and 5 for excitatory weights, and Ϫ5 to 0 for inhibitory weights. The final value of the weights was somewhat dependent on the initial conditions; and in general, the steady-state weights showed greater variability when starting from random initial conditions. However, the stability of the weights was unaffected. The other weights within the same segment (a, c, e, and f ) were kept fixed as in Table 1. All parameters were dimensionless, and time in the model was measured in arbitrary time units (t.u.).
Developmental time. For both activity-dependent models, simulations of development were conducted for a total time of 2 ϫ 10 5 t.u. The total simulation time depends on the time constant of weight update, the nature of spontaneous activity, and how spontaneous activity changes during development. The notion of "developmental time" in the model is useful in describing the gradual improvement of coordination in the activity patterns of the developing network.
Generating spontaneous activity. As there is little quantitative characterization of spontaneous activity in the Drosophila motor network during development, here we made simple approximations to generate spontaneous activity. Inputs of strength R i for duration T i were applied to each excitatory population E i (see Fig. 1). When the spontaneous input into a given population E i was sufficiently strong, E i became suprathreshold. At the time of threshold crossing, activity-dependent mechanisms become active to modify the weights. The duration of each spontaneous input, T i , was a random number drawn from a uniform distribution, U(x 1 , x 2 ), determining the total time for which the input of strength R i was applied to E i before choosing a new spontaneous input. Sampling the duration of spontaneous input from a distribution ensured that neuronal activity in different populations crossed threshold at different times as we did not want to introduce correlations in spontaneous activity. Because network connectivity is recurrent, the network itself generates a significant proportion of activity when the weights are sufficiently strong. This automatically generates correlations between the activations of neighboring segments, even though spontaneous activity itself is uncorrelated. The strength of spontaneous input, R i , was also a random number; this input was applied to E i for time T i . We used a truncated Gaussian distribution with mean and variance 2 , ͓N͑, ͔͒ ϩ , to ensure that spontaneous input remained positive.
The parameters of the distribution for T i did not greatly affect the frequency of suprathreshold events in a given excitatory population (data not shown). In contrast, R i had a strong effect upon the frequency of threshold crossings in the model network. This was partially because weights were modified precisely at the time when population activity crossed threshold. Therefore, we explored the effect of the distribution of R i on the outcome of weight development in the two plasticity models.
During development of Drosophila embryos, frequent random activity uncoordinated across segments gradually becomes replaced by coordinated activity ranging across several neighboring segments, with all contraction activity ceasing shortly before hatching when the larva emerges from the egg shell and starts to crawl (Crisp et al., 2008). This is consistent with a model where spontaneous activity decreases during development.
To model this effect, we considered decreasing either the mean or the variance of the Gaussian distribution for R i . In the Hebbian model, we decreased the SD in steps; however, the model was insensitive to how spontaneous activity was reduced. The homeostatic model was also robust to the decrease of spontaneous activity, although the speed of decrease of spontaneous activity affected the speed at which weights stabilized.
The Hebbian bidirectional model. Weight updates followed the ratebased formulation of Hebbian plasticity: if presynaptic precedes postsynaptic activity, then the corresponding weight is potentiated or else the weight is depressed (see Fig. 2 A, B). For any two populations connected by a weight, the activity that crossed threshold first was presynaptic. For instance, if E iϩ1 crossed threshold before E i , then b i f was potentiated, whereas b i b was depressed. To implement potentiation, a small constant, ⌬ ϩ Ͼ 0, was added to the weight as follows: and similarly, if the constant for weight depression was denoted by Unless otherwise noted, and b i f weakened: For the inhibitory weights, if I i crossed threshold before E iϩ1 , then d i f ͑t͒ strengthened: whereas if E iϩ1 crossed threshold before I i , then d i f ͑t͒ weakened: And, similarly, if I iϩ1 crossed threshold before E i , then d i b ͑t͒ strengthened: whereas if E i crossed threshold before I iϩ1 , then d i b ͑t͒ weakened: For a weight to be modified, activity of the presynaptic and postsynaptic populations connected by the weight had to be simultaneously suprathreshold. Weight changes were applied immediately after the second of the two populations connected by the weight crossed threshold. However, the model was robust to this requirement, and the weight could also be updated as the activity of each population dropped below threshold. The Hebbian efficacy model. In the Hebbian efficacy model, synchronous threshold crossing of activity in neighboring populations resulted only in weight potentiation and no depression. To ensure that weights did not grow without a bound, weight increase was limited by the frequency of contraction of the neuronal populations connected by the synaptic weight. Such activity-dependent control of the weights was motivated by the regulation of activity-dependent de-pression in developing motor networks (Fedirchuk et al., 1999;Crisp et al., 2008).
To implement this activity-dependent synaptic depression, a "synaptic efficacy" variable was introduced in the Hebbian model (Froemke and Dan, 2002). This efficacy variable depressed to 0 at the time when the activity of a neuronal population crossed threshold (see Fig. 2C). Then, in Equations 4 and 6, when the excitatory activity in two populations crossed threshold, the excitatory weights, b i f and b i b , were modified by ⑀ i f ͑t͒⌬ ϩ and ⑀ i b ͑t͒⌬ ϩ , respectively. Similarly, instead of Equations 8 and 10, the inhibitory weights, d i f and d i b , were strengthened by i f ͑t͒⌬ ϩ and i b ͑t͒⌬ ϩ , respectively, with efficacies as follows: Here, t e denotes the time when activity of the two neuronal populations connected by the weight crossed threshold (see Fig. 2C). The efficacy was held at 0 for a refractory time period t ref , after which it recovered to 1 with a time constant of rec . The refractory period was implemented to ensure that, after the weights had strengthened sufficiently to generate waves, threshold crossings of the activity in neuronal populations during propagating waves did not further modify the weights. The value of rec was chosen such that the period between two consecutive waves was larger than the refractory period t ref .
To equalize the frequency of threshold crossing among the different excitatory populations, an additional "end drive" was applied to the end excitatory populations, E 8 and E 1 (see Fig. 5). This drive was applied to one of E 8 or E 1 chosen randomly for a duration T i ϭ 3 t.u. instead of the spontaneous input R i . We assumed that the strength of this end drive increased as spontaneous activity decreased during development, reaching a final value of 1.7.
The homeostatic model. The homeostatic model was inspired by the Bienenstock-Cooper-Munro theory of synaptic plasticity (Bienenstock et al., 1982). Weight modification was bidirectional (potentiation or depression) to maintain a constant averaged activity of each excitatory population, r i , as follows: where r i ͑t͒ ¡ r i ͑t͒ ϩ 1 when E i ͑t͒ crosses threshold. (13) The time constant r was chosen to be much longer than the timescale of neuronal dynamics. Weight potentiation or depression depended on whether the average excitatory activity in a given segment was greater or smaller than a threshold for synaptic modification, . The homeostatic model is activity-dependent; thus, in the absence of spontaneous activity ( ¡ 0), the weights in the network will not be modified. The threshold is defined to be a function of the averaged excitatory activity itself: where r 0 ϭ 2.0 was a constant and p Ͼ 1 (here, p ϭ 1.5). The modification function depends on r i (t) and i ͑r i ͑t͒͒ as follows: Then, synaptic weight change was implemented according to the following: for the forward excitatory weights, and for the backward inhibitory weights for which for the backward excitatory weights, and for the forward inhibitory weights for which E iϩ1 is postsynaptic (i ϭ {1, 2,…, 7}). The constant was a small positive number to ensure that weights change slowly relative to the timescale of neuronal dynamics (e.g., 0.001). In the limit of r i ͑t͒ Ͻ Ͻ r 0 , This leads to ͑r i ͑t͒, ͑r i ͑t͒͒͒ Ͼ 0 because p Ͼ 1, and potentiation of the weights. In the limit of r i ͑t͒ Ͼ Ͼ r 0 , This leads to ͑r i ͑t͒, ͑r i ͑t͒͒͒ Ͻ 0 and depression of the weights. In both models, the weights were prevented from changing sign in agreement with Dale's law (Strata and Harvey, 1999).

Results
Can spontaneous activity successfully tune a bidirectionally connected recurrent network with weak connectivity by using local activity-dependent mechanisms that modify connection strength to produce functional output? Here we studied the generation of unidirectional (forward and backward) propagating waves of activity resembling those observed during crawling in the Drosophila larval motor network. Two main principles provided the basis for our developmental model: (1) the model should use local activity-dependent mechanisms to tune connection strength; and (2) the connections should achieve a stable configuration at the end of development such that the network generates propagating waves of activity with regular timing relationships (interburst intervals and duty cycles) in different segments. We used a recurrent network with bidirectional connectivity, which can generate forward and backward propagating waves when appropriately tuned ( Fig. 1) (see also Gjorgjieva et al., 2013;their Fig. 3). The model was initialized with weak connectivity between different segments incapable of generating waves. Weak connectivity in the network is likely specified by activityindependent mechanisms, which guide axons toward their correct target partners. These guidance mechanisms are widespread throughout the central nervous system (Dickson and Gilestro, 2006;Huberman et al., 2008;Klein and Kania, 2014;Hand and Kolodkin, 2015). Each excitatory population in the model network was triggered by patterned spontaneous input (see Materials and Methods), and depending on its strength and duration, increased the activity of the populations in random segments above threshold at different times. We sought to determine the nature of activity-dependent plasticity rules that modify the strength of excitatory and inhibitory weights, b and d, connecting neighboring segments using the spatiotemporal activation patterns of different segments. Having as a goal to produce robust propagating waves of activity with precise timing relationships, we compared a Hebbian with an activity-dependent homeostatic model for connection tuning based on local activity patterns.

The Hebbian bidirectional model does not generate functional weight distributions
Weight modification in the Hebbian model was based on Hebb's mechanism (Hebb, 1949): if the population activity in two neighboring segments is simultaneously above threshold, the weights between the neuronal populations should increase ( Fig. 2 A, B). Even though spontaneous activity in the model was uncorrelated across segments, correlations between the activations of neighboring segments exist because of the network's recurrent connectivity. Hebb's idea has been translated into various forms of correlation-based rules for synaptic modification, which commonly assume that synapses change strength in proportion to the correlation of presynaptic and postsynaptic activity. This includes spike timing-dependent plasticity, where the magnitude of synaptic potentiation and depression depends on the order of presynaptic and postsynaptic action potentials, as well as their timing (Markram et al., 1997). Hebbian plasticity has been successfully applied to activity-driven refinement of developing circuits and learning in neuronal networks in many brain regions, including hippocampus, neocortex, and cerebellum (Bliss and Lomø, 1973;Bliss and Collingridge, 1993;Malenka and Nicoll, 1999;Martin et al., 2000;Malenka and Bear, 2004). Correlation-based Hebbian rules have been often used to model map formation in the visual system (Miller et al., 1989;Miller, 1992). Therefore, we first characterize the role of Hebbian rules in connectivity tuning in the recurrent networks for motor output studied here.
A common challenge that Hebbian mechanisms face is unbounded growth of synaptic weights. Although upper bounds on synaptic strength can be imposed, and most biological systems probably have a saturation constraint limiting synaptic growth, it is unlikely that this upper bound is always fixed to the same value and that the system knows this value a priori. To prevent unbounded synaptic growth, we considered two implementations of the Hebbian model. One implementation is based on bidirectional plasticity, which evokes synaptic potentiation if presynaptic activity comes before postsynaptic activity, and depression if the order of activity is reversed (as in spike timing-dependent plasticity). This model has been referred to as the Hebbian bidirectional model (see Materials and Methods). The Hebbian bidirectional model fails to produce stable bidirectional connections because the weights develop in a biased way: some potentiate, whereas others decrease to zero, uncoupling neighboring populations (Fig. 3).
What is the reason for this bias? Because random spontaneous input is applied to each excitatory population in the network with initially small b and d, some weights randomly potentiate more than others. However, each pair of neighboring excitatory populations is connected by weights for wave propagation in both directions, so whenever the weight for one direction potentiates, the corresponding weight in the other direction depresses. This creates a feedback loop, in which stronger weights produce more Dashed lines indicate that these connections are plastic during development. The weights in the forward and backward direction are identified by superscripts "f" and "b," respectively. Each segment receives spontaneous random input of strength R i for a duration of T i , as illustrated on the left. B, An example of suprathreshold excitatory activity (threshold for detection of suprathreshold activity is 0.3) of each segment upon driving E 8 with P ext ϭ 1.7 for a network with parameters as in Table 1. Time is measured in arbitrary time units (t.u.). Inset, Excitatory activity (blue) and inhibitory activity (red) for segment A8 with a threshold indicated by the dashed line.  frequent suprathreshold activity, further amplifying the same sets of weights, whereas the weights in the opposite direction decay to zero (Fig. 3). Therefore, it is not possible to strengthen synaptic connections in both directions between neighboring excitatory populations appropriately for wave generation because of the asymmetric integration window of the Hebbian model. This is consistent with various models of spike timing-dependent plasticity (Abbott and Nelson, 2000;Song and Abbott, 2001). Nonlinear plasticity rules can generate functional bidirectional connectivity, but such models are based on detailed spiking and voltage trajectories of individual neurons (Clopath et al., 2010).

The Hebbian efficacy model requires balancing of activity across the network
We therefore considered an alternate implementation of the Hebbian model where synaptic potentiation is induced if presynaptic activity occurs before postsynaptic activity, augmented by activity-dependent synaptic efficacy to prevent unbounded weight increase (Fig. 2C), but without any synaptic depression. This model has been referred to as the Hebbian efficacy model.
The amplitude of weight modification in the Hebbian efficacy model is inversely proportional to the frequency of threshold crossings in segments connected by that weight (see Materials and Methods). The model produces bidirectional weights that become ordered in the direction of wave propagation (Fig. 4A). Thus, of all excitatory weights for forward propagation, the most posterior (b 7 f ) becomes the largest, whereas the anterior weights (b 1 f Ϫ b 3 f ) become the smallest. We found this to be due to a combination of (1) the difference in total input (synaptic and spontaneous) received by each excitatory population, and (2) the preference for synaptic potentiation of the forward weights induced by the order of excitatory segment activation, a hallmark of Hebb's rule. Indeed, segment A8 has only one neighbor and receives overall less drive than the other segments in the network; thus, its efficacy is the largest, evoking the highest potentiation of b 7 f . Although the most anterior segment, A1, also has only one neighbor, it can also be driven from its neighboring segment (A2) in the forward direction so b 1 f has a lower efficacy and does not potentiate as much. The weights at the anterior end and the middle of the network are most frequently activated; they have the lowest efficacy and are potentiated the least.
We examine this in greater detail. The threshold crossing frequency of excitatory activity (necessary to trigger weight change) depends strongly on the strength of the weights at a given time point during the development simulation (Fig. 5A). When weights are small, threshold crossing is purely the result of spontaneous drive, which is approximately uniform across all segments; thus, all segments exhibit approximately the same number of threshold crossings (Fig. 5A, blue). When the weights are sufficiently strong to generate waves, then the activation of any one segment propagates to all remaining segments, activating all of them uniformly. Hence, all segments again exhibit approximately the same number of threshold crossings (Fig. 5A, red). The biggest discrepancy between the frequency of suprathreshold activity in different segments occurs at intermediate synaptic strength (Fig. 5A, yellow). In this scenario, the end segments of the network, which have only one neighboring segment, receive the lowest synaptic input and so have the lowest overall input. Weights connecting neuronal populations whose activity crosses threshold more frequently due to the higher overall drive (in the middle of the network) have smaller efficacies compared with the weights at the ends of the network. Therefore, the former weights potentiate less than the latter. The outcome of this process produces variable steady-state weights (Fig. 4A). Taking the steady-state weights at the end of the developmental period fails to produce propagating waves of activity across the network (data not shown).
If Hebbian-style activity-dependent mechanisms indeed tune synaptic connectivity in the bidirectionally connected network, then to equalize the total amount of input received by each segment, the end segments must receive additional spontaneous input (Fig. 5B). Furthermore, as spontaneous activation of all segments decreases during development (Crisp et al., 2008), then the intensity of additional spontaneous drive to the end segments must increase. The increase in the end drive during development might correspond to the biological system progressively testing more frequently whether the network is fully assembled to produce waves. Thus, toward the end of development, the only drive to the network is to the end segments, which will initiate either forward or backward waves (Fig. 5B).

Weight distributions from the Hebbian efficacy model produce variable waves even with balanced activity
The Hebbian efficacy model with additional end drive successfully produces stable and functional weight distributions (Fig. 4B). To examine the progression of network output during development, we recorded the spatiotemporal activity patterns in the network over a given period (Fig. 6A-E) at the developmental time points denoted in Figure 4B (arrows). The network output gradually becomes more coordinated similarly to motor development in Drosophila (Crisp et al., 2008). Furthermore, partial motifs of coordinated output, which involve the simultaneous suprathreshold activity of excitatory populations in several neighboring segments, typically originate at one end of the network. This makes a concrete experimental prediction for how coordinated activity emerges gradually during development if the Hebbian efficacy model with additional end drive organizes the network.
The Hebbian efficacy model with balanced activity across the segments generated stable weight distributions, even when the initial conditions for the weights were randomly distributed in a given range, albeit still weak (Fig. 4C). Increasing that range of initial conditions fails to produce stable weights (data not shown). To determine whether the steady-state weights can indeed generate waves, forward waves were initiated by driving the most posterior population E 8 in 20 networks where the weights developed from random initial conditions (Fig. 6F ). The resulting waves have variable interburst intervals and variable duty cycles across the different segments that often fall outside the experimentally measured ranges (Fig. 6G,H ).

Wave sensitivity to spontaneous activity patterns in the Hebbian efficacy model
How sensitive is weight development and wave generation in the Hebbian efficacy model to changes in the properties of spontaneous activity? Modest changes to the distribution for generation of spontaneous input strength R i preserve the interburst intervals but retain the irregular duty cycles during wave propagation (Fig. 7). Some cases, however, produce a stable configuration of weights too small to generate robust waves (Fig. 7, bottom). More substantial variations in the spontaneous input strength result in more significant variability of the final weight distributions and ultimate failure of the network to generate waves. Similarly, stable weight distributions and propagating waves are produced only for a small subset of efficacy time constants and recovery parameters. Mistuning these parameters either destabilizes the weights so that they never saturate to a stable value or stabilizes the weights too early before they are sufficiently strong for wave generation (data not shown).

Conclusion for the Hebbian model(s)
For our networks with predetermined connectivity spanning only neighboring segments, we conclude that the basic requirement for coincident presynaptic and postsynaptic activity in the Hebbian models is insufficient to generate functional weight distributions that produce propagating waves with the appropriate timing relationships. Although model performance improved with the introduction of new assumptions (synaptic efficacy and additional spontaneous drive to the end segments), the final weight distributions and waves were highly variable. This suggests that Hebbian mechanisms likely play a minor role in tuning the motor network for wave propagation.

The homeostatic model
When experimentally challenged with perturbations in synaptic structure and function, neurons have the remarkable ability to regulate their synaptic strengths to bring activity back to a normal range. Various homeostatic mechanisms that maintain neuronal stability have been identified in invertebrates and vertebrates (Pérez-Otaño and Ehlers, 2005;Davis, 2006;Marder and Goaillard, 2006;Turrigiano, 2008;Pozo and Goda, 2010;Marder, 2012). One of the bestdefined systems for analysis of the homeostatic mechanisms that regulate synaptic efficacy is the Drosophila neuromuscular junction. The neuromuscular junction exhibits a strong homeostatic response to changes in postsynaptic excitability (Petersen et al., 1997;Davis et al., 1998). In contrast to the popular view that homeostatic plasticity is a slow phenomenon involving many neurons simultaneously, neurons may also undergo rapid synaptic tuning (Pozo and Goda, 2010). Therefore, we considered a different model to regulate synaptic strength: a homeostatic mechanism that maintains network activity at stable levels based on spontaneous input. In our implementation of the homeostatic model, synaptic weights are potentiated or depressed to maintain a constant average activity of the excitatory postsynaptic population in each segment (see Materials and Methods). Spontaneous input is necessary for synaptic change. As such, our homeostatic model falls under the category of local activitydependent models, which require that each neuronal population computes a running average of its activity and the weights are modified to maintain that quantity (Zenke et al., 2013). A different class of homeostatic models are nonlocal models, which require that all incoming weights to a neuronal population are kept normalized, even in the absence of activity (such as synaptic scaling) (Turrigiano et al., 1998).
The mechanisms for bidirectional homeostatic plasticity (simultaneously inducing potentiation and depression) have recently been demonstrated for the Drosophila neuromuscular junction (Gaviño et al., 2015), and also in Drosophila central neurons. The motor and visual system of developing larvae also shows bidirectional structural homeostasis (Tripodi et al., 2008;Yuan et al., 2011), as well as the mushroom body of the adult (Kremer et al., 2010). In our homeostatic plasticity model, weights potentiate or depress to maintain excitatory activity in all segments (postsynaptic to both excitatory and inhibitory weights) at the same target level. For instance, weights potentiate or depress based on whether the average excitatory activity r i in a segment i is larger than a modification threshold, ͑r i ͒; the threshold itself is a nonlinear function of the slow average of postsynaptic activity, as proposed in the Bienenstock-Cooper-Munro plasticity model (Bienenstock et al., 1982) (Fig. 8 A, B). To determine the steady-state activity, we can solve for the case when the weights are also at a steady state. Then from Equation 15, The steady-state excitatory activity is as follows: When r 0 ϭ 2.0 and p ϭ 1.5, then r eq ϭ 8.0 (Fig. 8A). Varying p affects the value of the steady-state excitatory activity r eq (Eq. 21) and the steady-state weights in the network (Fig. 8C). For instance, when p ϭ 1.7, maintaining r eq ϭ 5.4 requires smaller  ϭ 2). B, Same as in A but with additional end drive applied to one of the end segments chosen randomly with equal probability (reaching 1.7 at the end of the simulation). Spatiotemporal patterns of excitatory activity at the times denoted with arrows are shown in Figure 6. C, Same as in B but with random initial conditions starting in the range of 0 and 5. In all cases, spontaneous activity was generated with R i ϳ [N(0.5, 0.8)] ϩ and T i ϳ U(2, 3). Every 8000 time units, the SD ( ϭ 0.8) of the distribution for R i was decreased by 0.05. Only excitatory weights in the forward direction are shown as per Figure 3.
weights in the network (Fig. 8D). Similarly, when p ϭ 1.3, maintaining r eq ϭ 20.2 requires larger weights (Fig. 8E). The timescale r also has an effect on the steady-state weights. Increasing r results in lower average excitatory activity, r i (Eq. 13), and smaller steady-state weights (Fig. 8F ). Therefore, the parameters in the homeostatic plasticity model determine the average activity of the excitatory populations and the strength of the steady-state weights in the network necessary to maintain that activity.    Figure 4B. A, Scale bar applies to A-E. F, Applying P ext ϭ 1.7 to E 8 in the network with the final set of weights generates a forward wave. G, H, Duty cycle and interburst interval (mean Ϯ SD) for waves that were generated over 20 trials of wave development as in Figure  4C.   Figure 8 shows that, even at a constant level of spontaneous network activity during a simulation of development, the weights stabilize to maintain the target level of postsynaptic activity. If the amount of spontaneous activity decreases during the developmental period as the network assembles (Crisp et al., 2008), then the weights increase more gradually toward their steady-state value (Fig. 9A), with a time constant combining the slow synaptic weight change and the decrease in spontaneous activity. Weight evolution goes through a plateau before the network has fully assembled. This plateau occurs as a result of the interaction of spontaneous input and recurrent activity in the network when the weights strengthen sufficiently to propagate activity to neighboring segments. Changing the properties of spontaneous activity can change when (and whether at all) this plateau occurs (see Fig. 11).

The homeostatic model produces stable weight distributions
During forward wave propagation, all but one of the excitatory weights (b 2 f , …, b 7 f ) have as postsynaptic activity the excitatory activity of a segment in the middle of the body (E 2 , …, E 7 ). The postsynaptic activity for b 1 f is E 1 , which is coupled only to a single neighboring segment, E 2 . Thus, to achieve the same target activity, b 1 f potentiates more strongly than the other forward excitatory weights (Fig. 9A). The homeostatic model generates stable and functional weight distributions, even when the weights start randomly distributed in a given range (Fig. 9B,D). In contrast to the Hebbian model, the homeostatic model successfully generates stable weight distributions for larger ranges of initial conditions. For especially large ranges of initial weights, functional connectivity does not emerge as some weights continue to potentiate, whereas others depress, depending on their initial strength (data not shown).
Last, the homeostatic model can maintain weights at a functional range, even when network connectivity is perturbed. Following appropriate weight development in our functional model network, we "cut" all connections between segments A3 and A4 (by setting the weights to 0), which eliminated wave generation (Fig. 9E). Allowing the homeostatic model to act in the presence of spontaneous input recovered functional weights and appropriately timed waves (Fig.  9C,F).

The homeostatic model generates appropriately timed propagating waves
To examine the output of the model during development, we recorded the spatiotemporal patterns of activation in the network ( Fig. 10A-E) at the developmental time points denoted in Figure  9A (arrows). Coordinated output gradually improves in the network similar to the gradual improvement of motor output during Drosophila development (Crisp et al., 2008).
Are the stable weight configurations produced by this model able to generate robust propagating waves with appropriate timing relationships? Probing 20 networks, where the weights were generated from random initial conditions, shows that the homeostatic model can indeed generate waves with regular interburst intervals and duty cycles across the different segments that closely match experimental variability (Fig. 10G,H). This demonstrates that the homeostatic model successfully tunes network connectivity to a functional state, which can generate propagating unidirectional waves with the appropriate segmental timing relationships.

Wave sensitivity to spontaneous activity patterns in the homeostatic model
The amount of spontaneous activity, determined by the properties of the distribution of input strength into each excitatory population, determines the steady-state value at which the network weights saturate. To further investigate how the steady-state weights depend on the nature of spontaneous activity, we examined weight development and wave generation as a function of the mean () of the distribution for spontaneous input strength, R i ϳ [N (, )] (similar results were obtained when varying the SD; data not shown).
Increasing raises the average spontaneous input strength, R i , thus increasing the excitatory activity in each segment, E i . Because E i receives two types of input (spontaneous and synaptic input), to maintain an average activity at the target level r eq , the synaptic input from neighboring populations must decrease; this can be achieved by decreasing the steady-state strength of the weights (Fig. 11). Waves generated by the network using the final steady-state weights become longer as increases, and excitatory activity remains above threshold for a shorter time in each segment, resulting in shorter duty cycles and longer interburst intervals. Increasing Ն 0.5 makes the influence of spontaneous input on excitatory activity so strong that very little drive from neighboring segments is needed to maintain the average excitatory activity at target. Therefore, the weights become too small for the network to generate waves (data not shown). and consequently the steady-state feedforward excitatory weights. E, F, Varying r also affects the steady-state weights, even when the target activity is the same. Here, spontaneous activity was generated with R i ϳ [N(0.5, 0.4)] ϩ and T i ϳ U(2, 3). Activity remained constant during development. C-F, Weights shown as per Figure 3, but only forward excitatory weights for illustration.
This suggests that, for the homeostatic model to produce stable and functional network weights, spontaneous activity must operate within a range. It should be sufficiently frequent to drive excitatory activity above threshold and to enable weight potentiation; yet it should not be too frequent to allow steadystate population activity to be maintained by synaptic input from neighboring populations. The model is robust to changes in the amount by which spontaneous activity decreases during development, and it can also generate stable and functional weights even when spontaneous activity is constant during development (Fig. 8).

Conclusion for the homeostatic model
These results demonstrate that, based on the properties of spontaneous activity, different stable weight configurations can be achieved with the homeostatic model by maintaining a target level of postsynaptic activity for each weight. The weights can be reliably reproduced over multiple simulation runs and initial conditions, and no additional assumptions are necessary to bound synaptic growth as for the Hebbian model. A network with the final set of weights can generate propagating waves with precise timing properties, such as regular interburst intervals and duty cycles across the segments. This suggests that differences in the nature of spontaneous activity during development may underlie variability of network connectivity and the resulting wave properties produced by the final network. Therefore, developmental differences in spontaneous activity (endogenous or environmental) represent one plausible way by which such variability is achieved in the motor network of Drosophila (Berni et al., 2012;Gjorgjieva et al., 2013).

Discussion
How are neural circuits organized to maintain stable function and produce robust output? This task is made especially difficult during circuit development, when the properties of different circuit components are immature and changing, and input from the environment unreliable. Here we studied a model network of recurrently connected excitatory and inhibitory neuronal populations segmentally repeated in a single dimension, motivated by the motor network for Drosophila larval crawling. A key feature of this network, when appropriately tuned, is its ability to generate robust forward and backward waves of activity by driving each end of the network (Gjorgjieva et al., 2013). Examining sets of connectivity parameters that enable wave generation in the network has revealed nontrivial solution regions. We sought to identify plausible mechanisms by which these parameters emerge.

Activity-dependent tuning of network connectivity
Although there is ample evidence that spontaneous activity patterns instruct refinement of circuit connectivity in sensory systems (Huberman et al., 2008), the importance of spontaneous activity for the development of motor circuits is conflicting (Marder and Rehm, 2005;Myers et al., 2005;Li et al., 2007;Roberts et al., 2014). Network models of spontaneous activity in developing motor circuits have typically addressed the generation and properties of high-activity episodes interrupted by quiet periods in the embryonic chick spinal cord (Tabak et al., , 2001Marchetti et al., 2005;. However, unlike sensory systems, it is unclear whether this activity helps refine developing motor circuits. Imaging of spontaneous activity during embryogenesis in Drosophila larvae has revealed a gradual progression of motor output, and manipulations of this activity point to the role of activity-dependent mechanisms in the tuning of the network (Crisp et al., 2008(Crisp et al., , 2011Giachello and Baines, 2015). Using computational modeling, we have demonstrated that spontaneous activity can indeed tune a weakly connected recurrent network if the appropriate activity-dependent tuning rules are used. We examined two styles of activity-dependent mechanisms for weight development that through gradual improvement in network output lead to stable weight distributions to generate propagating waves with regular interburst intervals and duty cycles.
Hebbian-style models modify synaptic strength based on coincident presynaptic and postsynaptic activity. Although Hebbian mechanisms instruct activity-driven refinements in the developing visual system (Huberman et al., 2008), applying similar mechanisms for the tuning of a motor network for wave propagation yields less successful results. Unlike sensory systems where neurons topographically project from an input to a target layer, the motor network in Drosophila larvae is segmentally organized and recurrently connected. Although correlated activity emerges between neighboring segments, even when spontaneous input is uncorrelated, there is no obvious notion of presynaptic and postsynaptic: the same segment can be both presynaptic and postsynaptic to different synaptic connections. Thus, Hebbian mechanisms are less likely to play a fundamental role in fine-tuning network connectivity. Indeed, both Hebbian (bidirectional and efficacy) models fail to produce func-  Figure 10. B, Same as in A but with random initial conditions starting in the range of 0 and 5. C, Starting with the final weights in B (T 0 ), only those between segments A4 and A3 were cut (weights set to 0) as indicated by the arrow (T cut ), and allowed to develop again under the homeostatic model (T end ). D, Forward wave generated using the final weights in B at T 0 . E, The network fails to generate waves when the connections are cut using the initial weights in C at T cut . F, Forward wave generated by the recovered network using the final weights in C at T end . In all cases, spontaneous activity was generated with R i ϳ [N(0.3, 0.8)] ϩ and T i ϳ U(2, 3). Every 8000 time units, the SD ( ϭ 0.8) of the distribution for R i was decreased by 0.04. Only excitatory weights in the forward direction are shown in A-C as per Figure 3.
tional bidirectional weights that generate waves with appropriate timing relationships.
In neural systems as diverse as the crustacean stomatogastric system and the vertebrate visual system, homeostatic mechanisms control neuronal function through the regulation of synaptic efficacy and the modulation of intrinsic ionic conductances (Pérez-Otaño and Ehlers, 2005;Davis, 2006;Marder and Goaillard, 2006;Turrigiano, 2008;Pozo and Goda, 2010;Marder, 2012). Several forms of homeostatic plasticity have been proposed to stabilize network activity, such as regulation of the strength of synaptic transmission and synaptic scaling as a function of global network activity (Turrigiano et al., 1998;Turrigiano, 1999;Turrigiano and Nelson, 2004;Gonzalez-Islas and Wenner, 2006). Altered activity patterns can scale synaptic connections through homeostatic mechanisms also in motor circuits (Borodinsky et al., 2004). Our homeostatic model was based on the goal to maintain a target level of excitatory postsynaptic activity, which was achieved by modifying the synaptic weights for which the target population is postsynaptic. Spontaneous activity was necessary to drive weight modifications (Zenke et al., 2013).
Without any additional assumptions, the homeostatic model achieves stable and functional bidirectional weight distributions that generate waves with regular duty cycles and interburst intervals.
Moreover, different wave properties can be produced from different spontaneous activity patterns, providing a way to generate wave variability in the experimental system. Therefore, our results suggest that activity-dependent homeostatic mechanisms are more likely to tune weak bidirectional connectivity in a recurrent network than Hebbian mechanisms. During tuning, the network gradually improves the coordination of its output, matching the gradual emergence of coordinated output in developing Drosophila larvae. The resulting network can produce appropriately timed network-wide waves of activity that propagate in either direction despite bidirectional connectivity. This property can be matched to the ability of Drosophila larvae to crawl forward and backward with similar properties (Gjorgjieva et al., 2013). The homeostatic model can also restore functional weights following perturbations of connectivity in the presence of spontaneous input with similar patterns as used for developmental tuning. Thus, if activity-dependent homeostasis continues to operate after development to keep patterned activity robust, then the network must continue to receive external input. This is plausible because the Drosophila larval motor network receives descending input during continuous crawling, as well as additional proprioceptive input from the environment.  Figure 9A. A, Scale bar applies to A-E. F, Applying P ext ϭ 1.7 to E 8 in the network with the final set of weights generates a forward wave. G, H, Duty cycle and interburst interval (mean Ϯ SD) for waves that were generated over 20 trials of wave development as in Figure 9B. The model is compared with experimental data as in Gjorgjieva et al. (2013).