Terminal Schwann cell and vacant site mediated synapse elimination at developing neuromuscular junctions

Synapses undergo transition from polyinnervation by multiple axons to single innervation a few weeks after birth. Synaptic activity of axons and interaxonal competition are thought to drive this developmental synapse elimination and tested as key parameters in quantitative models for further understanding. Recent studies of muscle synapses (endplates) show that there are also terminal Schwann cells (tSCs), glial cells associated with motor neurons and their functions, and vacant sites (or vacancies) devoid of tSCs and axons proposing tSCs as key effectors of synapse elimination. However, there is no quantitative model that considers roles of tSCs including vacancies. Here we develop a stochastic model of tSC and vacancy mediated synapse elimination. It employs their areas on individual endplates quantified by electron microscopy-based analyses assuming that vacancies form randomly and are taken over by adjacent axons or tSCs. The model reliably reproduced synapse elimination whereas equal or random probability models, similar to classical interaxonal competition models, did not. Furthermore, the model showed that synapse elimination is accelerated by enhanced synaptic activity of one axon and also by increased areas of vacancies and tSCs suggesting that the areas are important structural correlates of the rate of synapse elimination.

. Developing sternomastoid NMJs at P0, P3, P7, and P16. (a,c,e,g) Montaged 2D electron micrographs spanning the entire NMJs in individual thin sections of mouse sternomastoid muscle fibers at P0, P3, P7, and P16, respectively. The nerve terminals of the NMJs contain synaptic vesicles and come within ~50 nm of the surface of a muscle fiber (M). The endplate on the muscle fiber is occupied by tSCs (green dotted line), vacancies (black dotted line), and axons (dotted line in color different from green and black). Note that the postsynaptic membranes at P0 and P3 are smooth but show multiple indentations called junctional folds (JFs) at P7 and P16. Multiple axons at the NMJs are also associated with tSCs; vacancies are present, some of which extend to the surfaces of the muscle fibers. (b,d,f,h) Segmentation of the muscle fiber surfaces (red), terminals of axons (here in shades of blue), tSCs (green), and vacancies at the same NMJs. Red line represents the muscle fiber surface. Scale bar, 1 µm.
P3, and P16 obtained by serial-section electron microscopy from a previous study 7 were used while those of 5 NMJs in the muscles at P7 were newly obtained by serial block-face scanning electron microscopy and examined for this study. The images in Fig. 1a,c,e,g show that each endplate is occupied by multiple axons, tSCs, and vacancies. Axon nerve terminals (purple, cyan and blue in Fig. 1) are outlined in purple, and tSCs surrounding the axons are present at all endplates (color-coded in green in Fig. 1b,d,f,h). Vacant sites (vacancies) devoid of contacts of either tSCs or axons are marked by black lines (Fig. 1a). The postsynaptic membranes of the muscle fibers are outlined in red and exhibit junctional folds (JFs) by P7, a feature of normal development and mature NMJs.
The three-dimensional (3D) surface models of the tSCs (green), vacancies (black), and axons (other colors) were generated after they were traced or segmented (Fig. 2a-h) as in the previous study 7 . The sites of the tSCs, axons, and vacancies opposed to their endplate were identified (See Fig. 1a) and their contact areas were measured by computing the entire area of each contact site using their surface models generated after segmentation (See Methods). The measured areas were normalized by dividing the areas by the total areas of tSCs, vacancies, and axons to compare their relative contact areas. From P0 to P3, the relative contact area of tSCs in the endplate increased on average whereas that of axons decreased and the relative area of vacancies changed little (Fig. 2i) as seen in the previous study 7 . We also found that the relative contact area of tSCs changed little between P0 and P3, and then decreased (Fig. 2i). The dominance of tSC occupation from P3 to P7 is consistent with its active role in synapse elimination by competing against axons to occupy the territory of the endplate proposed in the previous study 7 . In contrast, the relative contact area of vacancies was found to markedly decrease from P3 to P7 and then change little by P16, whereas the relative contact area of axons significantly increased from P3 to P16 (Fig. 2i). The increase in the axonal area suggests that axons pursue territory of the endplate more actively than tSCs and vacancies. Considering the involvement of Schwann cells and axons in the competition 7,26 , these dynamic changes are expected to be closely related to the process of synapse elimination. However, even with current advanced imaging technologies, direct observation of these dynamic changes is extremely challenging to quantify. Accordingly, we set out to take a modeling approach to investigate the relationship of the dynamic change in the relative areas of tSCs, vacancies, and axons within the endplate during synapse elimination.

A model of synapse elimination involving direct roles of tScs and vacancies. Previous modeling
studies have shown the transition of a polyinnervated synapse to a singly innervated synapse is a result of competition among different axons in an endplate 8,36 . However, the models did not address the presence of tSCs and vacancies at endplates 7 . To address their roles, we simplified the contact sites of tSCs, vacancies, and axons in each endplate so that the transition of the endplate is represented as the reorganization of their simplified contact sites (Fig. 3a), and based on the previous studies 7, 26 , we assumed that tSCs, vacancies, and axons compete against each other (Fig. 3b). Then a synaptic site formed by an axon has three different transition probabilities: a probability of transitioning into a vacancy (P AV ), a probability of transitioning into a tSC (P AS ), and a probability of no transition (P AA = 1-P AV -P AS ). Similarly, a synaptic contact site formed by a tSC has three different transition probabilities: a probability of transitioning into a vacancy (P SV ), an axon (P SA ), and no transition (P SS = 1−P SV −P SA ). A vacant site also has three different transition probabilities: a probability of transitioning into a tSC (P VS ), an axon (P VA ), and no transition (P VV = 1−P VS −P VA ). Therefore, there are six independent transition probabilities (P AV , P AS , P SA , P SV , P VA , and P VS ); however, measuring such transition probabilities experimentally with sufficient observations is currently impossible. So we took a computational modeling approach here. We first constructed the initial layout of the endplate composed of multiple contact sites made by tSCs (green), vacancies (black), and axons (other colors) using the empirical data of synaptic footprints gathered from serial images of mouse sternomastoid muscles at P0 (See Figs. 1 and 2). These initial layouts, as seen in Fig. 4a,f,k, depict the relative contact areas of the sites of tSCs, vacancies, and axons measured in the P0 data sets and the random distribution of those sites. As noted above, it is reasonable to assume that those sites compete to take over other sites in the endplate. However, unlike axonal competition models of synapse elimination 8 if those individual sites interact with each other with random transition probabilities, synapse elimination would proceed unreliably because it cannot account for the maintained presence of tSCs and vacancies after synapse elimination is complete. To test this possibility, we generated an idealized endplate randomly occupied by multiple contact sites of tSCs, vacancies, and axons. We kept the area ratios nearly those measured at P0 (See Methods) and assumed that the values of P AV , P AS , P SA , P SV , P VA , and P VS are randomly assigned from 0 to 1 so that a randomly chosen site can be transitioned into a site of a tSC, a vacancy, or an axon adjacent to the chosen site. The competition process was iterated by repeating the random selection of a transition site until there is no change in the composition of the sites. The simulation beginning with these initial layouts using random transition probabilities showed that the competition can lead to an endplate entirely occupied by only one type of axon after losing all other types of sites at the endplate sharply; however, the competition did not show the stable copresence of tSCs, vacancies, and a single axon after synapse elimination (Fig. 4). Moreover, other simulations led to an endplate filled only with tSCs or a completely unoccupied endplate (S1). Even with nonrandom, equal transition probabilities, synapse elimination did not reliably produce an endplate filled with one axon, tSCs, and vacancies where the values of P AV , P AS , P SA , P SV , P VA , and P VS were all one-thirds (S2 and S3). Although some degrees of randomness and/or equal degrees in the interactions might be present during synapse elimination, our simulation results demonstrated that such randomness or equal fitness is insufficient to reliably recapitulate synapse elimination. Accordingly, we assumed that tSCs, vacancies, and axons have different fixed transition probabilities at least during a portion of this period of time. To determine their constant transition probabilities, we developed a stochastic model of synapse elimination that employs the areas of tSCs, vacancies, and axons at different postnatal days and Markov Chain methods that are widely used in biological modeling studies 29,[31][32][33][34][35] . Our model assumes that tSCs and axons compete against each other as proposed from a previous study 7 , but here we newly propose that the competition of tSCs and axons is mediated by vacancies, and different types of axons compete to take over their adjacent vacancy in an endplate (Fig. 5). Specifically, the model assumes that vacancies mediate the transition between axonal and tSC sites without their direct transition (Fig. 5),  The regions of the axons, tSCs, and vacancies in proximity to their muscle fiber or footprints at their endplate were segmented and using the series of their segmented lines the contact regions of the axons, tSCs, and vacancies were marked in color using the same color-code on their muscle fiber (red). (i) The average relative areas at P0, P3, P7, and P16. From P0 to P3, the relative contact area of tSCs in the endplate increased from 31.0% to 57.0% on average whereas that of axons decreased from 51.6% to 24.8% on average and also that the relative area of vacancies changed little (17.4% at P0 and 18.3% at P3 on average). The relative contact area of tSCs changed little from 57.0% at P3 to 55.6% at P7 and then decreased to 40.7% at P16. The relative contact area of vacancies markedly decreased from 18.3% at P3 to 4.55% at P7 and then changed little (5.82% at P16) whereas the relative contact area of axons significantly increased from 24.8% at P3 to 39.8% at P7 and then to 53.5% at P16. www.nature.com/scientificreports www.nature.com/scientificreports/ and different axon types compete to occupy nearby sites. These assumptions allow us to describe the stochastic competition process as a Markov process 29,34 . Accordingly, we have the following equations: AA AV VA VS There are only 3 independent transition probabilities (P AV , P SV , and P VS ) because P AA , P SS , and P VA depend on P AV , P SV , and P VS based on the assumed relationship in Eqs. 1-3. In addition, it is assumed that each synaptic site can be occupied by a neighboring axon consistent with other and our studies 7,8,26 . Thus, when an individual axon occupies multiple neighboring sites near a recently vacated site, the axon has a higher probability of reclaiming that vacancy than other axons with fewer sites proximal to the vacated site. As in the previous section, we constructed the initial layout of the endplate composed of multiple contact sites made by tSCs (green), vacancies (black), and axons (other colors) using the empirical data at P0. Again these initial layouts (Fig. 6a,f,k) depict the relative ratio of the contact sites' areas measured in the P0 data sets and randomly distribute tSC, vacancy, and axon sites in an endplate. Then, the competition process among them was iterated by repeating the Figure 3. Schematic diagram of competing tSCs, vacancies, and axons. (a) An endplate covered with four different types of axon terminals (light blue, purple, cyan, and dark blue) and a tSC (green) transitions into an endplate covered with three of the axon terminals after removal of one axon terminal from the four and a vacancy (black) that takes the place of the removed axon. (b) Competition of tSCs, vacancies, and axons. A blue circle represents a synaptic site formed on a muscle fiber by an axon. A green circle represents a synaptic site formed on a muscle fiber by a terminal Schwann cell (tSC). A black circle represents a vacancy having no axon or tSC on a muscle fiber. All the three different synaptic sites can transition from one kind of synaptic sites to another. A contact site formed by an axon has three different transition probabilities: a probability of transitioning into a vacancy (P AV ), a probability of transitioning into a tSC (P AS ), and a probability of no transition (P AA = 1-P AV -P AS ). Similarly, a contact site formed by a tSC has three different transition probabilities: a probability of transitioning into a vacancy (P SV ), a probability of transitioning into an axon (P SA ), and a probability of no transition (P SS = 1−P SA −P SV ). A vacant site also has three different transition probabilities: a probability of transitioning into a tSC (P VS ), a probability of transitioning into an axon (P VA ), and a probability of no transition (P VV = 1−P VS −P VA ).  . Schematic diagram of vacancy mediated competition between tSCs and axons. A model of competition between axonal and tSC sites assumes that tSC-axon competition to occupy the territory of the endplate is mediated by their adjacent vacancies. Accordingly, a synaptic site formed by an axon has two different transition probabilities: the probability of transitioning from an axonal site into a vacancy (P AV ) and a probability of no transition (P AA = 1−P AV ). Similarly, a synaptic contact site formed by a tSC has two different transition probabilities: a probability of transitioning into a vacancy (P SV ) and a probability of no transition (P SS =1−P SV ); a vacant site that are not an axonal site either a tSC site has two different transition probabilities: a probability of transitioning into an axon (P VA ) and a probability of transitioning into a tSC (P VS = 1−P VA ). (2019) 9:18594 | https://doi.org/10.1038/s41598-019-55017-w www.nature.com/scientificreports www.nature.com/scientificreports/ random selection of a transition site in the layout until only one axon (identified by color) remained. As shown in Fig. 6, the model with fixed transition probabilities, determined using the measured relative areas of axons, tSCs, and vacancies at P3, reliably reproduced synapse elimination. The results showed that the number of different axons decreases exponentially consistent with other studies 3,8,15 . The results also showed that the contact sites of tSCs, vacancies, and axons in an endplate interchange dynamically as different axons compete against each other (Fig. 6); when synapse elimination is complete, the endplate is occupied by tSCs, vacancies, and only one axon ( Fig. 6i-k). The rate of synapse elimination is known to vary between different muscles and even between fibers in the same muscle 14 . Consistently, when we repeated the simulation 100 times, the number of iterations varied broadly (the inset of Fig. 6k); the average iteration to complete synapse elimination was 10300 ± 6719 (SD). Furthermore, direct measurements of the relative ratios of tSCs, vacancies, and axons at P3 were different from one endplate to another (Fig. 7a). Accordingly, it is reasonable to expect that our model using the different ratios may exhibit variable rates of synapse elimination. As expected, simulation results of our model with the different ratios of the 8 measured endplates showed that synapse elimination is complete at very different numbers of iterations, thus variable rates of synapse elimination (Fig. 7b). Despite such broad variation, we noted that the least number of iterations to complete synapse elimination is less on average with large relative areas of tSCs and/or vacancies than those with small relative areas (Fig. 7b) suggesting that tSC and vacancy occupation are closely related to the rate of synapse elimination.
Synapse elimination correlated with relative areas of tSCs, vacancies, and axons. To examine the correlation of the rate of synapse elimination with the relative areas of tSCs, vacancies, and axons with regard to their interdependence, we combined the relative areas into a composite area by dividing the product of the area of tSCs and that of vacancies by that of axons. The composite area was found to be negatively correlated with the least number of iterations to complete synapse elimination (Spearman Correlation, p < 0.05) as shown in Fig. 7c; the www.nature.com/scientificreports www.nature.com/scientificreports/ same correlation analysis using manually adjusted relative areas also showed that the composite area is negatively correlated with the least number of iterations to complete synapse elimination (Spearman Correlation, p < 0.05) as shown in S4 indicating that the composite area is positively correlated with the rate of synapse elimination. It is also reasonable to expect that the rate of synapse elimination decreases over time because the synapse elimination is complete in several weeks after birth. As expected, the simulation using the average ratio at P7 showed that synapse elimination is complete after 15700 iterations on average, which is markedly greater than that at P3 (S5). The average least numbers of iterations to complete synapse elimination between P3 and P7 were found to be significantly different (n = 100, t-test p < 0.05) indicating that the rate of synapse elimination is slower at P7 than P3. Similarly, simulations using average ratios at P16 showed completed synapse elimination after 15000 iterations on average, also greater than that at P3 (S6). The average least number of iterations to complete synapse elimination were significantly different between P3 and P16 (n = 100, t-test p < 0.05) indicating that the rate is less at P16 than that at P3. Our results are consistent with our expectation of a decreasing rate of synapse elimination over time. The average iterations of the complete synapse elimination between P7 and P16 were not significantly different (t-test, p = 0.59), which may suggest that the rate has already reached its minimum around P7.
Synaptic activity dependent selective and accelerated synapse elimination. It is widely held that synaptic activity is closely related to normal synapse elimination; several experimental studies support that axons having higher synaptic activity tend to remain at their endplate [37][38][39][40] . To test whether our model can predict this activity dependent synapse elimination, we assumed that among 9 different axons, one, selected at random, is more active than the others and the active axon sites have two-fold lower probability of having a site that undergoes the transition than other sites (See Methods). We generated the initial layouts using the average relative ratios of tSCs, vacancies, and axons at P0 and determined their transition probabilities using the average relative areas at P3 in the same way with the previous sections. Our simulation results showed that the active axon was the winner with the probability of 97%, consistent with the studies on synapse elimination and synaptic activity. Furthermore, the results showed that synaptic activity accelerates the rate of synapse elimination four-fold compared to that of control NMJs (Fig. 7d). When two axons were set to be more active than the others, our simulations predicted that the winning axon came from the active axons with the probability of 98%. Additionally, two active axons further accelerated synapse elimination but only by 10%, indicating that multiple coherently active axons in an endplate weakly promotes competition leading to little or no additional acceleration of synapse elimination (Fig. 7d). Increasing the number of active axons slowed synapse elimination and the rate became even less www.nature.com/scientificreports www.nature.com/scientificreports/ than that when no axons were more active (Fig. 7d). This demonstrates that our model can account for synaptic activity dependent axon selection, and the acceleration and delay of synapse elimination it predicts are consistent with other experimental studies 37,39,[41][42][43] .

Discussion
This study developed a stochastic model of synapse elimination that, for the first time, incorporates classical interaxonal competition with the recently proposed competition between tSCs and axons, predicting that the competition is mediated by vacancies, a third endplate occupant. With serial electron micrographs from 25 different developing mouse NMJs, the model employed the relative contact areas of tSCs, vacancies, and axons on muscle fibers as key parameters to determine the degrees of their interaction. Simulation results of our model reliably exhibited coexistence of axons with tSCs and vacancies during synapse elimination and even after synapse elimination was complete. Furthermore, our model successfully simulated expedited synapse elimination with synaptic activity, and predicts that the rate of synapse elimination varies during development suggesting that the ratios of the contact areas of these structural components are important to understand the roles of the components in synapse elimination, the mechanisms underlying impaired and enhanced synapse elimination, and the development of neural systems.

interaxonal competition and competition between axons and tScs for their adjacent vacancies.
A competitive process among axons innervating the same muscle fiber is thought to play an important role in synapse elimination. At P0 each axon occupies a similar proportion of the endplate area, so it is impossible to determine which axon will remain after synapse elimination by comparing their areas at this stage. Within days, however, some axon sites are vacated and the vacant sites taken over by different axons leading to the eventual loss of all but one axon. Accordingly, it is reasonable to assume that a vacancy can form randomly by the removal of an axon from the endplate of a muscle fiber. Any axons in proximity are expected to compete for the newly vacated site. Such competitive processes have been observed in experimental studies and examined theoretically in several modeling studies 8,15,22 . Specifically, in a model based on evolutionary graph theory, synaptic rearrangement occurs in a neuromuscular junction filled with several different axons; after random removal of sites, neighboring axons with equal fitness compete for the vacated sites, the model successfully simulated the competitive process of synapse elimination 8 . However, a role for Schwann cells, which occupy the endplate alongside axons, in synapse elimination 7,26 is not considered in the model. Since the model is designed to allow only the winning axon to remain after synapse elimination, the coexistence of tSCs and vacancies with it is neither accounted for nor possible (See Figs. 5, S1, S2, and S3). While previous studies provided indirect evidence for Schwann cell-axon competition during synapse elimination and propose a qualitative model of Schwann cell mediated synapse elimination 7 , this and other studies did not provide a quantitative model to characterize such a direct role for Schwann cells. With respect to vacancies as mediators of competition, we assumed random occurrence of a vacancy in a muscle fiber endplate similar to random competition models of different axons for synapse elimination 8 . Furthermore, we extend the random competition to tSCs by assuming that a vacancy can form by removing an axon or a tSC from an endplate. It is reasonable to assume that a vacancy can be formed by removing a tSC from an endplate with the concomitant competition for the vacated site if a random competitive process occurs between tSCs and axons 7,26 . We assumed that their transition probabilities are closely related to their relative contact areas as determined using their relative areas based on Markov Chain methods, increasingly used to understand the probabilistic behavior of the dynamic stochastic mechanisms associated with biological processes 34 . According to our model, axons and Schwann cells do not have equal fitness for taking over vacancies, and assignment of equal or random fitness rendered synapse elimination unreliable (See S1, S2, and S3). A recent study reported that the relative area occupied by tSCs at P3 is greater than the area at P0 and also that at P16 7 , proposing a mechanism of synapse elimination in which tSCs directly compete against axons to take over an endplate. A study by others examining NMJs of soleus muscles of P7-P8 mouse reported that tSCs have different Ca 2+ responses to weak and strong axons, and proposed another mechanism whereby tSCs differentiate weak axons from strong axons and are involved in selectively removing weak ones, supporting the active role of tSCs in synapse elimination 44 . Furthermore, a recent study reported that overexpression of NRG1-III enhanced tSC activity and expedited synapse elimination while reduced NRG1-III slowed synapse elimination 26 . Those studies raise the possibility that increases in the relative area of tSCs may speed up synapse elimination. When we modified the relative area of tSCs in our models, simulations reflected an enhanced rate of synapse elimination, consistent with the study. Furthermore, our model predicted that an increase in the relative area of vacancies would enhance synapse elimination more significantly than increasing tSC area alone. Based on the results, it is expected that relative areas of vacancies are large during the most active phase of synapse elimination. As expected, our measurement of the areas showed that the large relative area of vacancies maintained until P3 whereas the large relative area of tSCs appears between P3 and P7. Such differences may be closely related to the level of tSC activity and tSC-involved mechanisms of regulating the rate of synaptic elimination. As the relative areas of tSCs and vacancies change during synapse elimination, our model suggests that the rate of synapse elimination also varies accompanying the change.
Synaptic activity dependent synapse elimination. Synaptic activity is known to affect the rate of synapse elimination, and a number of studies have concluded that although synaptic activity itself is not a sole decisive factor of synapse elimination [45][46][47] , relatively active axons have competitive advantage 37,41,[48][49][50][51][52][53][54][55][56][57][58][59][60][61] . Presynaptic block by tetrodotoxin and botulinum toxin and postsynaptic block by α-bungarotoxin and curare led to the prolonging of synapse elimination [52][53][54][55][56] . Electrical stimulation of the nerves of 6-7 day old rat pups at 8 Hz for 4-6 hours per day for 2-4 days was found to decrease the level of polyinnervation from ~85% to ~50% 37 ; furthermore, electrical stimulation of muscles also was found to accelerate synapse elimination 39 . From these studies, it (2019) 9:18594 | https://doi.org/10.1038/s41598-019-55017-w www.nature.com/scientificreports www.nature.com/scientificreports/ is reasonable to expect that enhancing synaptic activity accelerates synapse elimination and provides the active axon with an advantage. Consistent with prior studies, our simulation results showed that an active axon is more likely to win sole occupancy of the endplate and that such synaptic activity accelerates synapse elimination. Furthermore, our model predicted that when two different axons have the same degree of high synaptic activity compared to other axons, the effect of accelerating synapse elimination is significantly weakened. Our model also predicted that when more than two different axon types have the same degree of high synaptic activity, synapse elimination slows rather than accelerates, suggesting that multiple coherently active axons delay synapse elimination instead of accelerating it. Although our simulated high synaptic activity of an axon or more is difficult to be converted into a specific stimulation frequency currently, our results are consistent with several other studies reporting that asynchronous synaptic activity promotes synapse elimination whereas synchronous synaptic activity delays it 41,46,57,58,60,61 .
As with any model of a biological process, it is impossible to describe all factors in play. For example, intrinsic differences in fitness among multiple axons may be present. Kasthuri and Lichtman proposed that motor neurons innervating fewer muscle fibers than others are more competitive than others 62 . Since our model predicts that the axons of the more competitive neurons tend to possess greater contact areas than others, measurement of the areas at a specific postnatal day such as P0 can be used to test the proposed idea by examining the distribution of the areas. The contribution of muscle fibers, such as neurotrophic factor release or removal of acetylcholine receptors 38,39 , may influence the process of synapse elimination. Increase or decrease of histocompatibility complex class 1 molecules and N-Methyl-D-aspartic acid (NMDA) receptors has been shown to accelerate or delay synapse elimination, respectively 63,64 . In the case the increase in their level is expected to enlarge the areas of vacancies and/or tSCs following the accelerated synapse elimination whereas the decrease in the level is expected to reduce those areas following delayed synapse elimination. But the intrinsic differences, contribution of muscle fibers, and other factors known to be implicated in synapse elimination are not considered in our model because there are no reliable structural correlates with them yet. Despite these limitations, our stochastic model (See Fig. 8) reliably simulates synapse elimination, agrees well with the observed coexistence of Schwann cells, vacancies and axons at endplates at its conclusion, and provides an explanation of the different dependences of synapse elimination on synchronous and asynchronous synaptic activities. Furthermore, the model reveals a novel role for vacant sites in synapse development closely related to the competition between tSCs and axons. Testing the predictions made by the stochastic model of synapse elimination at single and multiple endplates would contribute to our more complete understanding of synapse elimination and offer insights into the nature of synapse development and network refinement.

Methods
Tissue preparation and imaging for serial-section transmission electron microscopy. 20 different serial electron microscopy data sets from sternomastoid muscles of wild-type C57B/6 mice of either sex obtained from a previous study 7 were analyzed. The mice were obtained from the Jackson Laboratory, bred in a colony at the University of Texas at Austin and Texas A&M at College Station, and killed on the day of birth (P0; n = 10; 5 mice), the third postnatal (P3; n = 8; 3 mice), and the sixteenth postal natal days (P16; n = 2; 2 mice). Experimental procedures were approved by the Institutional Animal Care and Use Committees at the University of Texas and Texas A&M University and all experiments and methods were performed in accordance with the relevant guidelines and regulations. The generation of wild-type mice and their sample preparation for electron microscopy for mice of P0, P3, and P16 have been described previously 7,65,66 . Simply, animals were perfused transcardially with 0.1 M sodium cacodylate buffer, pH 7.4, followed by the same buffer containing 2% paraformaldehyde and 3% glutaraldehyde. Sternomastoid muscles were removed and fixed overnight at room temperature in the same fixative. Muscles were washed with cacodylate buffer and stained en bloc in 1% osmium tetroxide, 1% ferrocyanide in cacodylate buffer for 5 h, washed with water, and then stained in 1% aqueous uranyl acetate for 2 h. Muscles were dehydrated in graded alcohols and acetone, then embedded in Epon 812. Regions containing NMJs were located by thick (0.5-1 µm) sections made on a Leica Ultracut UTC Ultramicrotome with glass knives, stained with 1% toluidine blue, and examined under a light microscope. Blocks were further trimmed using a 45 degree cryotrim tool (Diatome) and serial 65 nm (silver) sections cut with a 35 degree diamond knife (Diatome). Sections were mounted on formvar-coated Synaptek slot grids (Electron Microscopy Sciences) and observed under a FEI Tecnai Spirit electron microscope with an AMT Advantage HR digital camera. NMJs were identified by the apposition of nerve terminals to within ~50 nm of the muscle fiber surface in the presence of Schwann cells. Digital images of individual NMJs and the immediately surrounding area were captured from each section at a magnification of 16,500×. At this magnification, ~6-12 overlapping images were acquired to include the entire extent of the NMJ. These individual images were then montaged by manual manipulation in Adobe Photoshop to create a single image. Such montaged images over 70-150 serial sections were then imported into the software program Reconstruct 67 and calibrated based on section thickness and known distances in the specimen. In software, the montage in each section was aligned linearly to each adjacent montage using 6-8 objects common in both and starting with the central montage in the stack and working toward each the first and last montage. Once aligned, structures of interest in each section were segmented (traced) using the tools provided in the software. Measurements and 3D renderings of these structures were then generated using Reconstruct.
Tissue preparation and imaging for serial block-face scanning electron microscopy. Tissues of a mouse at the seventh postnatal day (P7) were prepared differently and imaged using serial-block-face scanning electron microscopy (SBSEM). Tissue preparation and staining was done with compliance to Renovo Neural Inc. (Cleveland, OH). Briefly, P7 C57/Bl6 mice were transcardially perfused with 2.5% glutaraldehyde and 4% paraformaldehyde in a 0.1 M sodium cacodylate buffer. Sternomastoid muscles were then removed and stored in perfusion solution for 1 day prior to staining. Tissue was then stained with uranyl acetate and osmium-ferrocyanide, followed by graded-dehydration and embedded in Epon resin. SBSEM image data sets of postnatal day 7 NMJs  Fig. 1e. The tissue blocks were mounted, examined, and sectioned in a Zeiss Sigma VP scanning electron micro-scope equipped with a Gatan 3View in-chamber ultra-microtome stage with low-kV backscattered electron detectors optimized for 3View systems. NMJs were initially identified within cross-sections of muscle tissue by axon terminals in direct apposition to junctional folds on the muscle fiber membrane and capped by terminal Schwann cells. Regions of interest were chosen to include the NMJs in whole cross-section. The first sample block was sectioned cross-sectionally along Figure 8. A stochastic model of tSC and vacancy mediated synapse elimination. According to the model, a vacant site or vacancy in an endplate of a muscle fiber forms randomly by a removing tSC or nerve terminal of an axon. The transition probability of a tSC to a vacancy and an axon to a vacancy are determined by measured relative areas of tSCs, vacancies, and axons. Subsequently, the vacancy is taken over by adjacent axons or tSCs. But the probability of its taking over by a specific axon is positively correlated with the axon's relative contact area surrounding the vacancy. In other words, the greater the contact area of an axon surrounding the vacancy is, the greater the probability of its taking over the vacancy is. When the synaptic activity of an axon is greater than other axons in an endplate, the probability of the axon's taking over is greater than other axons leading to acceleration of synapse elimination consistent with previous and other studies. Furthermore, our model newly predicts that as the relative area of the vacancy or tSC in an endplate increases synapse elimination speeds up raising a testable hypothesis that both the measurable relative areas of vacancies and tSCs are important structural correlates of the rate of synapse elimination. (2019) 9:18594 | https://doi.org/10.1038/s41598-019-55017-w www.nature.com/scientificreports www.nature.com/scientificreports/ the long axis of the muscle fibers. A series of 465 EM images in a field size of 130 µm × 50 µm were acquired at 2.2 kV with a resolution of 8.1 nm per pixel and 90 nm per slice. Eight NMJs were captured at P7 with 5 of those NMJs fully reconstructed for synaptic contact areas.
A model of tSc and vacancy mediated synapse elimination and its simulation. We used circular contact sites assuming that all of the sites are composed of 9 different axons, which were the maximum number of different axons observed at P0, tSCs, and vacancies consistent with our observation of them from 2D scanning electron microscope images, and their initial layout is generated based on the recently reported ratios of contact areas of the tSCs, vacancies, and axons at P0. To simulate the initial layout, it is assumed that the total area is covered by tSCs about 30%, vacant sites about 16%, and axons about 54%, based on their recently measured ratios at P0. All of the sites are randomly formed within a round region of a diameter of 300 pixels, and each site has a diameter of 30 pixels as shown in Figs. 4 and 6 and all the supplemental figures. Vacant sites are black, tSC sites green, and axonal sites in different colors (blue, red, dark blue, yellow, pink, purple, dark green, gray, and brown). Thus, each of the 9 different axons takes up about 6% of the total area. Examples of randomly generated initial layouts at P0 are shown in Figs. 4a and 6a. We used the optimal probabilities accounting for the synapse elimination and all the ratios of different kinds of synaptic sites using Markov chain properties 30 assuming that the ratio of their contact areas with a muscle fiber in a NMJ reflects the efficiency of synapse elimination. We noted that at P3, the ratio of the contact areas of tSC sites greatly increased and gradually decreased later on raising a possibility that the ratio of the areas of tSC sites plays an important role on the process of synapse elimination consistent with the recent study 7 . Accordingly, the measured ratio of the contact areas of tSC sites, vacancies, and axonal sites at P3 was used here to obtain the optimal transition probabilities by determining the probabilities that can stabilize the ratio using Eq. 4 based on the Markov Chain properties 30 . where → r f , a vector having the final ratios. We have used the average ratios of the developing NMJs at different postnatal days ranging from P0 to P16 as shown in Table 1. Then, transition probabilities were adjusted to generate the final ratios close to those obtained at P3 from our previous study ( → r f 7 ;) using multiple NMJs of neonatal mice. According to Eq. 4, two independent variables of all the three independent variables can be expressed in terms of the remaining independent variable. Thus, there is only one undetermined transition probability in the matrix. We chose the transition probability from a vacancy to a tSC site (P VS ) as the undetermined probability. As P VS is varied, the rest of the transition probabilities are determined using the value of P VS and the ratios at P3. Then, for each specific P VS , the process of the synaptic elimination was iterated until it was complete. But it was noted that for some values of P VS a few other transition probabilities happened to be negative or undetermined. Thus, the value of P VS was set as 0.6, which allowed us to use all of our measured ratios of tSCs, vacancies, and axons such as those of multiple mouse NMJs at P3, P7, and P16 and ensured that all the other transition probabilities were reliably computed from our measured ratios; the range of the probabilities was broad from ~0.4 to ~0.7. For the specific P VS , the process of the synaptic elimination was iterated until it was complete or the number of iterations reached a specified maximum number of iterations that was sufficient for completing the synapse elimination (See Figs. 6-7, S5 and S6).
To simulate the synaptic activity dependent synapse elimination, the initial layout was constructed in the same way. By assuming that only a type of axons is active in an endplate, all the sites of the active axon among 9 different axons were selected and we set the probability of each of the sites to have a chance to undergo the transition to be two times less than each of all sites excluding the active axon; then as described above, the process of synapse elimination was iterated until it was complete or the number of iteration reached a specified maximum number of iteration that is sufficient for completing the synapse elimination (See Fig. 7). We also assumed that more than one type of axons are active in an endplate. In the case, we set the probability of having the transition site in each site of those active axons to be two times less than each of all sites excluding the active axons similar to the case of a single active axon (See Fig. 7d).
Simulations based on random transition probabilities. Because each transition probability might fluctuate randomly, first all of the probabilities are assumed to be random ranging from 0 to 1. In each step of the simulation, a contact site from all the sites in an endplate was randomly chosen. Then the site is converted into a different site with randomly assigned transition probabilities while the relationships of the probabilities as described in Fig. 3b and its legend are maintained. The step was repeated more than 10000 times to ensure that only one type of axonal site survived at most as shown in Figs. 4 and S1.  Table 1. Average ratios of the contact area of tSCs, mouse NMJs measured at P0, P3, P7, and P16 (n = 10, 8, 5, and 2, respectively).