A computational model of bidirectional axonal growth in micro-tissue engineered neuronal networks (micro-TENNs)

2.1Micro-TENN fabrication and experimental measurements

Micro-TENNs were generated as previously described [4]. Briefly, agarose (3% w/v) was cast in a custom-designed acrylic mold to yield microcolumns with an outer diameter of 345 or 398 μm and inner diameter of 180 μm. Microcolumns were UV-sterilized and cut to a specified length before the lumen was filled with an ECM comprised of rat tail collagen 1 (1 mg/mL) and mouse laminin (1 mg/mL) adjusted to a pH of 7.2–7.4 (Reagent Proteins, San Diego, CA). To create the neuronal aggregates, embryonic day 18 (E18) cortical neurons were isolated from rodents and dissociated. The resultant single-cell suspensions were added to custom PDMS pyramidal wells and centrifuged at 200 x g for 5 minutes to force the cells into spheroidal aggregates. Following 24 h incubation at 37°C/5% CO₂, aggregates were seeded within the microcolumns to generate unidirectional (with one aggregate) and bidirectional (with one aggregate at each end) micro-TENNs. Micro-TENNs were then grown at 37°C/5% CO₂ with half-media changes every 48 hours. To fluorescently label aggregates, adeno-associated virus 1 (AAV1) was sourced from the Penn Vector Core (Philadelphia, PA), packaged with the human synapsin 1 promoter and either green fluorescent protein (GFP) or the red fluorescent protein mCherry, and added to the pyramidal wells containing the aggregates (final titer: ∼3 ×1010). Aggregates were kept at 37°C, 5% CO₂ overnight before being seeded in micro-columns as described.

During the design and early development of the model, unidirectional and bidirectional micro-TENNs were generated with approximately 15-30E³ neurons per aggregate and lengths ranging from 2.0-9.0 mm (n = 39), with growth rates analyzed as described [5]. To identify aggregate-specific axons over time, a set of 3.0mm-long, bidirectional “dual-color” micro-TENNs were simultaneously generated such that one aggregate expressed green fluorescent protein (GFP) while the opposing aggregate expressed mCherry (n = 6). Finally, for quantitative validation of the growth model, 2.0mm-long, unidirectional micro-TENNs were transduced to express GFP and generated with approximately 20E³ neurons per aggregate (n = 6) or 8.0E³ neurons per aggregate (n = 6) for characterization as described below. Micro-TENNs were imaged under phase contrast microscopy (magnification: 10x) at 1, 3, 5, 8, and 10 days in vitro (DIV) using a Nikon Eclipse Ti-S microscope paired with a QIClick camera and NIS Elements BR 4.13.00 (National Instruments). In addition to phase contrast microscopy, the bidirectional dual-color micro-TENNs were imaged at 1, 2, 3, 5, and 7 DIV using a Nikon A1RSI Laser Scanning confocal microscope paired with NIS Elements AR 4.50.00.

To quantify micro-TENN growth rates over time, the longest identifiable axons were measured from phase images at each DIV using ImageJ (National Institutes of Health, MD). Lengths were measured from the leading edge of the source aggregate (identified at 1 DIV) to the neurite tip, and growth was measured until axons from the aggregate either spanned the micro-TENN length (unidirectional) or began to grow along axons from the opposing aggregate (bidirectional). Growth rates were averaged at each timepoint to obtain a growth profile for unidirectional micro-TENNs with 20E³ and 8.0E³ neurons/aggregate. The peak growth rates for each group were compared using an unpaired t-test, with p < 0.05 set as the baseline for statistical significance.

To characterize axonal density with respect to cell count, phase images of unidirectional micro-TENNs with either 20E³ (n = 6) or 8.0E³ (n = 6) neurons/aggregate at 5 DIV were imported into ImageJ. 10-μm long rectangular regions of interest (ROIs) spanning the inner diameter (final ROI dimensions: 180 μm x 10 μm) were taken at 50% and 75% of the micro-TENN lengths. The axon density at these two locations was quantified as the percentage of the ROI populated by axons. Densities were averaged for the 20E³ and 8.0E³ groups and compared at each location via unpaired t-test with p < 0.05 as the baseline for significance. All data presented as mean±s.e.m.

To characterize axon distribution, unidirectional micro-TENNs were fabricated and labeled with GFP (n = 5). At 10 DIV, micro-TENNs were gently drawn into a 22-gauge needle and vertically injected into a block of “brain phantom” agarose (0.6% w/v). Micro-TENNs were injected such that the aggregate was ventral with axon tracts projecting downward. Post-injection, micro-TENNs were imaged on a Nikon A1RMP+multiphoton confocal microscope paired with NIS Elements AR 4.60.00 and a 16x immersion objective. Micro-TENNs were imaged with a 960-nm laser, with sequential 1.2 μm-thick slices taken along the longitudinal axis (i.e. X-Y projections along the micro-TENN length). Post-imaging, the X-Y projections were used to generate a 3D reconstruction of the micro-TENN; cell bodies, axon bundles, and single axons were then manually identified via co-registration of the X-Y projections and 3D structure.

2.2Computational model development

This framework aims to emulate the growth and bifurcation of micro-TENN neurons by using simple diffusion principles. It avoids most of the underlying biological complexity (due to the lack of external guiding molecular cues) by considering only a few key parameters. Each tip of each neuron is considered as a diffusion source in free space. The elongation and the growth direction of each of the neurites in the model are guided by the contribution of the concentration gradients generated by of the rest of the neurite tips. The motivation for using a concentration gradient to mimic the growth behavior of neurons comes from our experimental observations and from various approaches in the literature. Mathematically, the concentration gradients help drive bundling, and from our experiments we see that axons bundle together as shown in Fig. 1. In addition, in the real brain, bundling of axonal fiber tracts is shown by diffusion tensor imaging [22]. The goal of micro-TENNs is to eventually be considered “replacement” fiber tracts that could be implanted. In both cases of real fiber tracts and micro-TENNs, bundling is observed suggesting an attractive mechanism that we can capture using a concentration gradient approach. Furthermore, it is known that axons have growth cones which are sensitive to various stimuli [23] that create attractive forces. Our model is a phenomenological construct of these ideas that run in a fast fashion ignoring the detailed molecular aspects but focusing on capturing the desired morphology. The bifurcation of the neurites is assumed to be a stochastic process, i.e., branching is associated with a time dependent probability function at each node. All of the tips of the neurite tree are assumed to participate in the extension and branching process; furthermore, extension and branching of each node are modeled as independent processes. This has computational advantages such as improved speed and ability to parallelize on a large scale.

The model uses a continuous space/discrete time approach to allow freedom in the outgrowth direction and elongation. Space is bounded by the inner diameter of the hydrogel micro-column, the diameter and length of the tubular hydrogels, 180 μm and 2 mm respectively, are based on experiments previously performed by the Cullen Lab [1, 3]. In the micro-TENNs, axonal extension is measured approximately every two days; as such, the size of the fixed time interval of the model is 1% of this two-day interval (i.e., 28.8 minutes). In each time step, each individual axonal tip may: (i) extend, (ii) bifurcate into two daughter branches and (iii) change growth direction. In the present implementation, the model uses fixed time steps with functions built upon the diffusion equation and concentration gradients for extension, turning, and branching. The model is developed with the following conditions: the extension rate and turning direction depend on the concentration gradients at the terminal segment of each axon, the extension rate decreases exponentially to zero value [12] as the neurites stop growing due to the limitation of space and essential biochemical factors [14], and branching probabilities grow as a function of the simulation time.

3.1Modeling setup: Diffusion equation and concentration gradient

Many stochastic models of neuronal activity are based on the theory of diffusion processes [24]. In our growth model, the tips of the neurons are modeled as diffusion sources in free space, assuming a constant isotropic diffusion coefficient:

The sum of all concentration gradients at each neurite tip is used to help guide the direction of growth for the next time step.

3.2Growth tip position

The coordinates at any given next step of a growth tip (P₂) are determined from the coordinates in the previous step (P₁), the outgrowth direction (D₂), and the extension rate (L); the new position in each time interval is given by:

3.3Direction of neurite outgrowth

The outgrowth of neurites is a complex process that is far from being fully understood. In actual biological processes, the outgrowth direction of neurites depends on many intracellular and extracellular cues, which may cause large fluctuations in outgrowth directions [22, 23]. Our model assumes that the new outgrowth direction depends on the previous outgrowth direction and on the concentration gradients of the growth tips. For each growth tip, the concentration gradients are normalized, and the new outgrowth direction is given by:

Fig.2

a) Visualization of Equation 4 which shows how the outgrowth of neurites is implemented. The variables are defined in the text. This controlling this term allows change in the magnitude of deviation of the growth direction. b) Pseudocode and schematic of the code implementation requires an initial parameter specification; parameters as sensitivities (S₁, S_d1), base growth rate (v₀), growth rate based on gradient strength (v_0grad), time constant for growth rate (tau), time constant for branching rate (tau_b), branching probability, bundling (RI), helix formation (chi), total steps, number of cells, inner radius of micro-column among others; these are specified based on experimental/known values. The algorithm then selects among four cases based on chi and bundle values, these cases are used to determine positions after the initial positions. Once calculated, the positions are plotted for its visual inspection and files with coordinates are saved for posterior visualization in Paraview.

3.4Rate of neurite extension

The rate of extension of a neurite (L) may vary considerably and is determined both by the external environment and by the internal state of the neurite in the micro-TENNs [24–31]. In general, the extension rate decreases gradually with increasing distance from the soma [15]; in our model, the description of neurite extension rate follows the trend of experimental growth rate measured in unidirectional micro-TENNs. In each time step, the elongation of a single neurite is represented by the function:

3.5Branching probability, rate of branching, and growth rate after branching

Neurite branching patterns are complex and show a large degree of variation in their shapes, random branching on the segment indeed results in large and characteristic variations in the structures of the tree; as previous research has highlighted [12, 32], branching is assumed to occur exclusively at terminal nodes. Our model describes branching as a stochastic process: for each time step, for each of the terminal nodes in the growing tree, a branching probability pb_j to form two new daughter nodes in a given time interval is assigned. The probability of a branching event at each given terminal node j is given by:

The time-dependent branching probability pb_j of a given terminal node j is dependent on several terms: the steady state branching probability (Pb|_t = ∞), simulation time t, the branching time step t_ib, and a branching time constant τ_b.

The latter equation assumes that the branching probability of terminal nodes at each time step remains constant for all tips but branching probabilities are growing with the total simulation time; such function was necessary to match the shape of increasing number of dendritic terminal nodes during outgrowth of the micro-TENNs. Whenever a branching event takes place, two daughter terminal nodes are instantaneously added to the end of the existing terminal segment [33], which then becomes an intermediate segment. The stochastic process of branching is also restricted by another random value E₃; branching could only occur when both pb_j and E₃ are greater than a certain value B, the value of B can be determined by the branching probability from experimental data.

The growth rate of the generated trees is closely related with segment outgrowth direction and extension and we only consider the extension distance in the axis direction as the growth distance; therefore, the growth rate is determined by the difference between the z components of the nodes. In the model, we force the growth cones to extend preferentially in the axis direction and the turning is relatively small, in this way the segment extension rate is the strongest factor to determine the growth rate. Estimates of v₀, and τ were obtained from the experimental growth rate and the experimental tuning of the elongation parameters which involves a comparison of the experimental and model segment extension rate.

As stated above, with the current implementation the branching probability increases with simulation time and the steady-state branching probability Pb and time constant τ_b are supposed to be extracted from experimental images. By controlling the values of Pb and τ_b, we have control over the morphology of the simulated neurites, since Pb controls the branching density and τ_b controls how early in the process branching begins, e.g., in the extreme case of P_b = 0 we can generate a morphology with no branching.

4.1Radial constraint

Micro-TENNs are grown within miniature tubular hydrogels experimentally, in the simulations the outgrowth process is also restricted within the tubular space (in this particular case, 180 μm in diameter [1, 3]); however, the model allows different simulation radii to be employed. At each time step, the radial components of all the terminal nodes are tracked, if a radial component of a given node does not satisfy the tubular constraint, the node will be re-oriented to stick on the tubular wall.

4.2Overlap avoidance

The branches and extensions of neighboring neurites often target a shared or adjacent position, i.e., many of the neurites are competing for space and avoiding overlap. Space competing is achieved as a result of the growth based on the concentration gradients explained above. At the same time, overlapping is avoided by checking the distance from the new positions to the surrounding existing segment tips; the model re-orients the growth direction whenever the extended position is sufficiently close to others.

4.3Bundling and helicity

One more feature included in the model is the capacity to accommodate fiber bundling. This is achieved by including an attraction term, i.e., the new position in each time interval becomes:

Helicity is another feature that was introduced to the model aiming to reproduce observed experimental micro-TENN morphologies. Very often, single axons and axonal bundles once reaching the inner wall of the micro-column, form a helix (Fig. 1); our model allows control over both the slope of the helix formed by an axon (or axonal bundle), as well as its helicity (handedness).

4.4Parameter tuning

Finding a best fit of the model-generated neuronal morphologies for an experimental data set requires an iterative comparison of experimental and model shape properties. When tuning the parameters of the model, some can be directly related to properties of the experimental data or images and can be obtained from them, e.g., the time step is selected to be 0.02 days, based on actual experimental extension measurements (see Results section). Similarly, parameters as v_0grad, v₀, and τ, which predict the growth rate in axis direction, are susceptible to the same strategy. Noticeably, since v₀ can be/is extracted from the experimental data, the selection of the value of the diffusion coefficient (D) is guided by the restriction:

Other parameters as τ_b, Pb and S₃, govern the branching and fully determine the topological structure of the generated trees, they are directly related to segment branching rate which can also be inferred/estimated from the experimental observations.

Axon growth and bundling along the tube

Slices of the model-generated morphologies were analyzed with ImageJ using the analyze tool after contrast filter/tuning to count to axons and measure the occupied areas. Individual axon areas were divided by the total area in each slice and the resulting ratios for each RI were analyzed as frequency distributions.

Implementation, scaling and computational performance

The code was implemented in Python 2.7-3, even though Python is not recognized as a fast computing language the runtimes are not excessive and the simulations can be performed on standard/modern computers. The algorithm (Fig. 2b) employed by our model is scalable; simulations were performed for 100, 1000 and 10000 cells for the unidirectional micro-TENN with diameter 180 μm and length 2 mm with runtimes were 1.99, 28.93 and 1500.45 seconds, respectively. The simulations were carried out on a laptop (Windows 10 Enterprise 64, Intel i7-7700HQ CPU @ 2.80 GHz, 2801 Mhz, 4 Cores, 16GB DDR3 RAM), all code with is freely available at: https://github.com/PSUCompBio/GrowthModel.

5.1Examples of micro-TENNs morphologies

The model simulates satisfactorily the growth processes for unidirectional and bidirectional micro-TENNs when grown to 2000 μm. Variation of model parameters like sensitivity to attraction (S₃) and radius of influence (RI) allows us to generate/reproduce different morphologies seen in experimental setups. Figure 3a shows various unidirectional morphologies modelled for micro-TENNs seeded with 3000, 6500 and 10000 cells; it demonstrates how different values of these parameters control the morphologies, e.g., S₃ takes continuous values between zero and one, with S₃ = 0 corresponding to no bundle formation and S₃ = 1 corresponding to tight bundle formation (Fig. 3b and 3a respectively). Higher RIs show fewer bundles for each cell density (Fig. 3a).

Fig.3

Example of unidirectional model-generated micro-TENN morphologies: (A) Twelve morphologies for micro-TENNs of length L=2000 μm and radius R=180 μm, number of seeded cells 3000, 6500 and 10000, and RI=70, 80, 90 and 100 μm, respectively. S₃ = 1 causes tight bundle formation. (B) S₃ = 0 guarantees no bundle formation.

The case of bidirectional morphologies is presented in Fig. 5, S₃ = 0 ensures no bundle formation in spite of the denser setup (Fig. 5) while intermediate values of S₃ show not so compact bundles (Fig. 4).

Fig.4

Example of unidirectional model-generated micro-TENNs: Micro-TENNs of length = 2000 μm and radius R = 180 μm, number of seeded cells 3000, 6500 and 10000, and RI = 90 μm. Here S₃ = 0.75, bundles are not compact.

Fig.5

Example of bidirectional model-generated micro-TENNs: Analogously to the unidirectional case, the model can generate bidirectional morphologies. The micro-TENNs have length L=2000 μm and radius R=180 μm, the number of seeded cells is 3000. Here S₃ = 1 and RI= 90 μm shows bidirectional axonal bundle formation (left), in the case S₃ = 0 no bundles are formed (right). Neuronal cell bodies (aggregate) are labeled on each end but are hidden for clarity.

5.2Experimental calibration of axonal growth rate

Growth rates were computed from the model and compared to those measured from experimental images (Fig. 6a,b) of unidirectional micro-TENNs of approximately 8000 and 20000 neurons; rates are not statistically different (Fig. 6c,d). The optimized parameters, v₀ = 15, v_0grad = 0.008 and E₂ (a random uniform value between 0.8 and 1) provided an excellent fit with the experimental data. The results show that extension rates increase up to three DIV and then decrease as the distance from the soma augments, the same trend can be seen for both the unidirectional and bidirectional growth rates (Fig. 6e,f).

Fig.6

Calibration of simulated growth rates in vitro. (A) Phase images of a unidirectional micro-TENN with approximately 8,000 neurons at 1, 3, and 5 DIV. (B) Phase images of a unidirectional micro-TENN with approximately 20,000 neurons at 1, 3, and 5 DIV. Both micro-TENN groups exhibited rapid axonal growth over the first few DIV. Scale bars: 100 μm. (C) Growth rates from both micro-TENN groups at 1, 3, and 5 DIV. Micro-TENNs with ∼20,000 neurons/aggregate exhibited qualitatively faster growth rates than those with ∼8,000 neurons/aggregate, although there were no statistically significant differences in growth rates. (D) Axon density at 5 DIV across the two groups at 50% and 75% along the micro-TENN length (as illustrated in (B)), quantified as the percentage of the microcolumn channel occupied by axons. Micro-TENNs with ∼20,000 neurons/aggregate showed higher axon densities than those with ∼8,000 neurons/aggregate, although this was only significant at 50% along the micro-TENN length (*** = p < 0.001). (E) Experimental versus model growth rate for unidirectional micro-TENN with 8,000 neurons. (F) Experimental versus model growth rate for unidirectional micro-TENN with 20,000 neurons.

Figure 7a shows a 3D-reconstruction of a unidirectional micro-TENN at 10 DIV. Four corresponding slices were extracted for cross-sectional comparison to the computed results. Figure 7b–e show each slice marked with arrows to distinguish between neuronal cell bodies, axon bundles, and single axons. In order to compare morphologies, we use two model-generated X-Y projections along the Z-axis (Fig. 8) that resulted from a micro-TENN simulation (inner radius 180 μm, length of 2 mm, number of seeded cells 20,000). The model generates realistic morphology with single axons and axon bundles along the inner wall comparable to the experimentally reconstructed morphology in Fig. 7.

Fig.7

(A) 3D reconstruction of a unidirectional, GFP-positive micro-TENN at 10 DIV. (B–E) X–Y projections of the micro-TENN from (A) at the sections outlined in red. Orientation of the z-axis (positive) is into the page. Neuronal cell bodies (arrows) can be seen near the aggregate region in (B), from which axonal bundles (triangles) project and split into individual axons (caret). Scale bars: 200 μm (A); 50 μm (B).

Fig.8

Model generated morphology of a unidirectional micro-TENN. (A) X–Y projection at Z = 910 μm; (B) X–Y projection at Z = 1380 μm. Inner diameter 180 μm, length 2 mm and approximately 20,000 neurons.

From different slices along the z-axis in our simulations we analyzed the percentage of area covered in the slice which serves as descriptor of the axon bundling along the tube. It should be noted that according to the models, the highest average size for the axons is reached between 660 μm and 990 μm in z-direction, i.e., where the bundles are bigger along the tube. Experimental confirmation is needed to test our observations. Different combinations of simulation parameters create particular micro-TENN morphologies, in an attempt to qualitatively characterize them we measured the area ratio of the bundle/micro-TENN cross sections at 5 equidistant z-values along the micro-TENN for 3 different RI (Fig. 9d,e,f). The resultant distribution for each RI is distinct and allows qualitative comparison between simulation morphology data and experimental microscopy micro-TENN data when it becomes available. The most abundant type are the smaller ratio areas (individual axons), both RI 50 and 70 have bigger counts of this type than RI 90; on the other hand, RI 90 presents the bigger area ratios (bundles). These results reflect the effect of RI in bundle formation using our model.

Fig.9

Bundle % area coverage (circles) and average bundle size (triangles) as a function of length in the micro-TENN. (A) RI = 50; (B) RI = 70; (C) RI = 90. Combined area ratio distributions from five z-cross sections from simulated microTENNs with (D) RI = 50, € RI = 70 and (F) RI = 90, respectively. The distinct distribution correspond to distinct microTENN morphologies.

6.1Future work

One major objective of further building/developing the Bidirectional Growth Model is to generate simulations of detailed neuron growth patterns to ultimately enable the study of functional connectivity that our research group has begun [36]. For instance, these neuronal growth patterns could be used as the input of a spiking model to study the firing patterns within micro-TENNs. The framework could also be used to extract detailed connectivity information in scenarios following implantation; patterns at the distal ends of micro-TENNs upon integration with the host brain neurons could be modeled/described/anticipated.

The neuronal growth patterns provide the information for searching locations of synaptic connections and help to establish the spiking network simulation. In biological neuronal networks, synapses form in regions where tissues are in sufficiently close proximity and synaptic connectivity is estimated (proximity criterion of 0.5 μm). According to experimental design, synapses occur close to the cells aggregates, which are in the 100 μm range from each end of the micro-TENNs. Figures 10A/2B give examples of the locations of these synaptic sites where high connectivity is presented. Further development of the model could introduce additional guidance and attraction/repulsion molecular cues, once such experimental information is available, thereby systematically adding complexity and the ability to capture synergistic and/or competing features of intrinsic and extrinsic growth parameters.

Fig.10

Growth pattern and synapses examples. Illustration of synaptic formation close to the cells aggregate. The spheres indicate the synapses location. We consider that a synapse is formed when 2 fibers are closer than a threshold distance of 0.5 μm. Synapses occur close to the aggregate (within 100 μm). The connectivity information is extracted for the spiking network simulation. (A) The output image of the model: synaptic formation in bidirectional micro-TENNs without axonal bundles. (B) The output image of the model: synaptic formation in bidirectional micro-TENNs with axonal bundles.