Overview
This review is largely about methods for modelling neuronal development, rather than the latest research results. We will, however, point out the major findings of various models along the way. The review is not exhaustive of all theoretical work on neuronal development (for a comprehensive overview see [13]), but concerns mathematical models aimed at understanding the biophysics of the outgrowth of neurites. Most models consider only specific aspects of neurite outgrowth, and so in turn we will look at models of (1) neurite initiation (2) neurite elongation (3) axon pathfinding, and (4) neurite branching and dendritic shape formation.
Neurite initiation and differentiation
The first problem concerns how neurites form from an initially spherical neuronal cell. Then, subsequent to initiation, there is a process of differentiation in which one of the neurites becomes the axon, and the remaining neurites form the dendrites. These stages of growth are illustrated in Figure 1. Mathematical modelling of these processes has been undertaken by Hentschel and colleagues over the last decade [14–18]. They investigated how initial neurite outgrowth could be established by a positive feedback mechanism triggered by small inhomogeneities in the cell surface [14, 15, 18]. Calcium is assumed to be the morphogen that promotes cell outgrowth. In an excitable membrane model, local membrane potential is modulated by the influx of sodium through voltage-dependent sodium channels. Calcium influx also occurrs through voltage-dependent channels. Small bulges in the cell surface produce focal depolarizations leading to elevated sodium and calcium entry and hightened calcium-dependent outgrowth of the bulge, eventually leading to neurite formation. Sodium, calcium and voltage gradients are established from the distal to proximal ends of the neurite, in accord with experimental evidence.
The model was simulated using a cellular automata approach in which the cell interior is divided into 1 μm2 blocks and the cell surface is described by a linked list of 1 μm segments of membrane. Each segment contains information about the local transmembrane potential, the conductance of sodium and calcium and their submembrane concentrations. The concentrations of sodium and calcium in each cytoplasmic block are calculated on the basis of a quasi-equilibrium assumption. Growth proceeds by the addition or subtraction of segments from the cell surface. Full details of the simulation procedure are given in [18]. Further outgrowth of neurites is assumed to be rate-limited by the availability of a particular, but unspecified chemical at each neurite tip [16–18]. This chemical is produced in the cell body and transported by diffusion and active transport to the end of each neurite emanating from the cell body. This scenario was modelled by a set of coupled ordinary differential equations (ODEs) describing the chemical concentration in the cell body, c0, and in each of n neurite terminals, c1 to c
n
, and the lengths of each neurite, l1 to l
n
. The equations are:
where T
i
is the transfer rate of chemical from the soma to a neurite tip, D is the diffusion constant, A is the cross-sectional area of a neurite, V0 is the volume of the soma, V
i
is the volume of a neurite tip, and F is a growth-dependent active transport rate. The elongation rate of a neurite is proportional (α) to the concentration in the tip, c
i
, and chemical is consumed in the tip at a rate proportional (G) to the elongation rate.
A "winner-take-all" dynamic instability can emerge if one neurite has a slightly larger initial growth rate. This leads to greater consumption of the chemical at this neurite tip and the subsequent capturing of more chemical by this neurite than by the others. This quickly leads to rapid growth by this one neurite, with all other neurites having only very slow growth. This is characteristic of the differential growth of the neurite that becomes the axon in real neurons. This model is able to reproduce the results of a number of axotomy experiments in which severing of the nascent axon at different lengths, relative to the other neurites, results in either the axon reestablishing itself or another neurite becoming the new axon [18]. Such competitive effects between neurites have also been investigated in similar ODE models of the growth of branched dendritic structures [19–21].
Neurite elongation
Following initiation, neurites elongate and branch, reaching lengths of tens to hundreds of micrometres for dendrites, to over a metre or more for long axons. Modelling has been used to explore just how neurites elongate and the biophysical constraints on the lengths that may be achieved. This work has centred around the dynamics of cytoskeleton construction, particularly of the microtubules that form the major supporting scaffold within neurites [19, 22–24]. Elongation is assumed to be a function of the amount of available free tubulin at the growing tip of a neurite, which is assembled onto the ends of microtubules, extending the internal cytoskeleton. The constraints to growth are the production rate of tubulin and its transport from its site of synthesis, typically assumed to be the cell body (but see Alvarez et al [25]), to the neurite tip by diffusion and active transport. An illustration of this model is given in Figure 2. The most mathematically sophisticated treatment of this problem is the continuum model developed by McLean and colleagues [23, 24] in which a PDE describes the tubulin concentration along the length of a single, unbranched neurite. This model essentially consists of the following four equations:
The free tubulin concentration at a point x in the neurite at time t is given by c(x, t). The length of the neurite at time t is given by l(t). Space x lies in the domain 0 to l. Tubulin moves by active transport (a) and diffusion (D), and degrades with rate g. Synthesis of tubulin in the cell body at rate ∈0c0 results in a flux of tubulin across the boundary at x = 0. Assembly of tubulin onto microtubules results in a flux of tubulin across the boundary at x = l (assembly flux ∈
l
and return flux ζ
l
) and a change in length (assembly rate r
g
and disassembly s
g
).
This model is a generalisation of previous ODE and algebraic treatments of this scenario [19, 22]. Its solution, both analytically and numerically, is complicated by the fact that it is of the moving boundary type ie the spatial domain changes over time due to changes in the length, l. A stable and accurate numerical solution has been proposed in which the spatial domain is discretized into a fixed number of N grid points and integration is carried out on a spatially-transformed domain in which x lies only in 0 to 1 [26]. Transforming back into the real length domain results in the grid points growing further apart as the length increases.
Steady-state analysis has revealed that growth can proceed in three different regimes, as determined by the relative proportions of construction, due to the synthesis and active transport of tubulin and its assembly onto microtubules, and dissipation, due to tubulin diffusion, degradation and its disassembly from microtubules [23, 24]. If construction dominates dissipation, then large growth results. If dissipation dominates construction then only short lengths are achieved. Moderate growth ensues if construction and dissipation are approximately balanced. The transition between growth regimes is highly nonlinear. Growth in the large and short regimes is quite damped and is determined by active transport in the large regime and diffusion in the short regime. Oscillations in length can occur in the moderate regime where growth is affected by both the active transport and diffusion of tubulin [26]. The ODE models described so far only prescribe concentrations at a small number of points, typically the cell body and neurite tips (but see work on compartmental models, described below). Numerical solution of this PDE model enables the visualisation of tubulin concentration over space and time, to any desired resolution, as the neurite grows [23].
Aspects of this model have been the subject of more detailed treatments. Active transport of tubulin is here assumed to proceed with a flux proportional to the local concentration. In reality, tubulin dimers undergo periods of active transport when they are attached to the molecular motors, interspersed with periods of free diffusion. A PDE reaction-diffusion-transport model that covers this scenario, and considers both uni- and bidirectional motors for the active transport of intracellular organelles has been investigated by Smith and Simmons [27].
Modelling has also been directed at examining the detailed process of microtubule assembly. Monte Carlo simulation of the reaction and diffusion of individual tubulin dimers near the tip of a microtubule has been used to investigate the outgrowth dynamics of a microtubule and the sources of variability that might underly the observed dynamic instability in which a microtubule may switch from elongation to retraction and vice versa [28]. A mean field approach has been used to determine the diffusion-limited rate of microtubule assembly, with the conclusion that observed assembly rates are not limited by the diffusion of free tubulin [29].
The growth cone at the tip of a neurite also plays a significant part in neurite elongation and is not explicitly included in the microtubule assembly models. However, microtubules extend into the growth cone and interactions between the growth cone and microtubules are significant for neurite outgrowth. Microtubules exhibit so-called dynamic instability, the random alternation between microtubule growth and shrinkage. This process plays an important role in the motility of growth cones and their finger-like protrusions, the filopodia. The volume of a growth cone is so small that growing or shrinking microtubules are expected to cause fluctuations in the concentration of free tubulin. Monte Carlo simulations have shown that these fluctuations have a significant effect on microtubule dynamic instability [30]. They shorten microtubule growth and shrinkage times and change their distributions from exponential to non-exponential, gamma-like, which, interestingly, allow optimal search behaviour by microtubules [31]. Thus, growth cone volume, via its influence on the fluctuations in the concentration of free tubulin, can affect dynamic instability and, consequently, growth cone motility.
The growth cone exerts tension on the trailing neurite which influences microtubule assembly rates. Precisely how tension affects assembly has been the subject of a number of modelling studies (reviewed in [32]). Physical considerations essentially lead to a description in which the rate of microtubule assembly is an inverse function of the elastic tension, F, on the microtubules due to stretching of the neurite:
where l is microtubule length, r
g
is the unmodified assembly rate and s
g
is the disassembly rate (assumed not to be a affected by tension). Tension from the growth cone can relieve this elastic tension (lowering F) and thus promote assembly and outgrowth. In the PDE model described above, all tensions are assumed to be constant throughout. A description that considers variation of tension over time needs the inclusion of a hetergeneous environment that is detected by the advancing growth cone. Such models have been formulated for studying axon guidance, as described below.
The growth cone itself contains a complex actin cytoskeleton. It has been proposed that this cytoskeleton, when supporting lamellipodia at the leading edge of the cone, undergoes a caterpillar like movement to propel the growth cone along [33]. A more detailed modelling study of actin cytoskeleton dynamics in lamellipodia supporting cell movement is reported in [34].
Axon pathfinding
A fundamental feature of the outgrowth of an axon is its ability to follow a path to an appropriate target. This requires the axon's growth cone to sense environmental cues and to be able to turn towards the desired direction, or away from an incorrect location. Modelling efforts have been directed at the growth cone's ability to sense external chemical gradients [35], its turning ability [32, 33, 36], and the formation of nerve tracts in which groups of related axons bundle together and travel to a remote target region [37, 38]. Related work, which will not be covered here, includes the problem of topographic map formation in which the topography of the sending region of neurons is preserved in the pattern of connections formed by the arriving axons over the target region [39].
The most fully worked example of axonal elongation and direction finding is that of Aeschlimann [32, 36].
She considers the outgrowth of a single, unbranched axon in a two-dimensional external environment containing chemical gradients and glial cells (Figure 3). Elongation and bending of the axon is determined by visco-elastic stretching and nonelastic extension due to microtubule assembly. Direction finding and steering is the function of filopodia extending from a circular growth cone at the axon tip. Filopodia elongate to a fixed maximum length, then retract. Retracting filopodia exert tension on the growth cone, which in turn influences axon elongation. The direction of elongation is determined by the vector of forces generated by the distribution of filopodia on the surface of the growth cone. The probability that a new filopodium will emerge at a point on the surface is determined by the local submembrane calcium concentration, which is influenced by existing filopodia. A filopodium can detect a chemical gradient along its length [35] or physical contact with an external glial cell, with both signals resulting in an influx of calcium into the growth cone at the base of the filopodium. This results in a positive feedback, with the calcium concentration being increased most around existing filopodia sensing the largest gradients, and in turn increasing the probability of new filopodia emerging in this area. Localised groups of filopodia contribute strongly to the tension vector and can thus cause the growth cone to turn towards and grow along the maximum external concentration gradient of an attracting chemical. This complex model contains many components and includes a coupled ODE description of neurite elongation and probabilistic descriptions of filopodium formation.
Maskery et al [40] use a much simpler model of axonal pathfinding in which a growth cone undergoes stochastic changes in direction or is repulsed by external cues to demonstrate that guidance by external cues is most effective when growth is at a transition between being dominated by stochastic or deterministic effects.
Models of how groups of growing axons form bundles, or fascicles, have been explored by Hentschel and van Ooyen [37]. In their most efficient (de)bundling model, a diffusible signal released by the target region provides long-range axon guidance. Bundling of the axons occurs due to a diffusible chemoattractant released by the axonal growth cones, causing the axons to grow towards each other. Debundling of the axons as they come near the target occurs due to a chemorepellant released by the axonal growth cones, with the amount of chemorepellant being a function of the concentration of the target-released attractant. Debundling must occur for the axons to spread out and innervate the full extent of the target region. The chemoattractant and repellant released by the axonal growth cones remain hypothetical. They could plausibly be mediated by contact interactions between growth cone filopodia, rather than be diffusible molecules [37].
The two dimensional implementation of this scenario involves treating each target cell and each growth cone as a point source for attractant or repellant. The quasi-steady-state concentration gradient of each chemical from each source is calculated and summed. The change in position, r
α
, of each growth cone, α, is then calculated as a function of the two attractive and one repulsive gradients:
where ∇ρ
x
is the concentration gradient of growth cone (c) or target (t) attractant, or growth cone repellant (r) and λ
x
is the associated growth rate constant.
Segev and Ben-Jacob [41, 42] have used a similar modelling approach to explore more general chemotactic effects in neuronal network formation in a 2D environment. In their models single axons grow out from each soma, with their direction of growth being determined by a variety of attractive and repulsive environmental cues. In particular, each soma is assumed to emit a repulsive signal. Once an axon reaches a specified length it emits a pulse of an attractive signal. A soma then responds by itself emitting an attractive signal if it detects an above-threshold concentration of growth cone attractor. As a consequence a growth cone is then attracted to the nearest soma. This model was used to investigate the patterns of network connectivity that could result from various initial spatial distributions of somas. In some simulations spatial distributions of glial cells that emitted chemotactic signals were also present. Validation of these models awaits new measures for characterizing and distinguishing between different network organisations.
Neurite branching and dendritic shape formation
Neurites do not just elongate, they also make repeated branches. Branch formation is either due to a bifurcation of the growth cone, or the interstitial formation of a new branch part way along an existing neurite. Growth cone bifurcation seems to largely underly the formation of many different types of dendrite [43]. Interstitial branching is important for the formation of axon collaterals [43] and also for the creation of short oblique branches in certain dendrites [44]. We will concentrate here on mathematical models that consider the biophysics of growth cone bifurcation and the implications for the creation of the dendritic structures that characterize different neuronal types. The models can be divided into those that describe external influences on the growth cone and those that consider internal constraints on cytoskeletal construction and stability. There is also a range of statistical models that are aimed at the reconstruction of realistic dendrites, without necessarily describing the developmental process (see for example [45, 46]). Such models will not be considered here.
Models of branching due to external influences consider tension on the growth cone due to filopodia sensing cues in different directions in the environment [47, 48]. If there is sufficient tension in opposing directions on different sides of the growth cone, then the cone splits apart to form two daughter branches, which proceed to elongate in the direction of the sensed cues (Figure 3). These models consider growth along a plane, equivalent to an isolated neuron growing in culture. The external environment is populated by attracting cues, which are either point sources or continuous chemical gradients. A number of filopodia radiate from a point growth cone, with direction and tension dependent on the pattern of external cues. A tension threshold determines when the growth cone branches as it elongates through the environment. The micro-details of the internal actin cytoskeleton, its stability and how it is rearranged following a branching event are not considered in these simple models. The output of these models is comparable to data on branching angles and segment lengths from cells grown in tissue culture [49].
Theoretical studies indicate that dendritic branching angles may follow a principle of neurite volume minimisation [50], or equivalently, equilibration of tension forces on growth cones [51]. Network formation may refine branch angles to minimise neurite lengths [52]. Statistical analysis also indicates a strong tendency for neurites to grow away from their parent cell body, thus influencing branching angles [53]. Models of internal constraints are similar to the models of elongation based on the construction of the microtubule cytoskeleton. In addition to considering the rate of neurite elongation due to microtubule assembly, these models also consider the stability of the resultant microtubule bundles as determining the likelihood for a neurite tip to bifurcate. As with the external models, an explicit model of the growth cone and its actin cytoskeleton currently are not included.
The basic model considers whether the production and transport of an unspecified branch-determining substance imparts constraints on branching [54]. Statistical models of dendrite formation indicate that the probability of a neurite bifurcating is modulated as the tree grows and is a function of the number of terminals in the tree and the number of branch points between a particular terminal and the cell body (centrifugal order) [3]. The basic model shows that modulation of branching probability as a function of the number of terminals and their centrifugal order can arise from variations in the availability of a branch-determining substance at each terminal due to the diffusion and active transport of that substance from its site of production in the cell body. Such a substance can plausibly be identified as tubulin, or possibly a microtubule-associated protein, such as MAP2, that influences the stability of microtubule bundles in dendrites.
More sophisticated models describe both elongation and the branching rate explicitly as functions of tubulin and MAP2 concentrations at a terminal (Figure 4) [55, 56]. Both tubulin and MAP2 are assumed to be synthesised in the cell body and transported by diffusion and active transport along the growing dendrites. At a terminal, the amount of free tubulin and MAP2, and the phosphorlyation state of MAP2 determine the elongation and branching rates. Dephosphorylated MAP2 acts as a cross-linker between microtubules, stabilising the bundles and promoting microtubule assembly. Phosphorylated MAP2, on the other hand, loses its cross-linking ability and destabilises the microtubule bundles, thus increasing the likelihood of a bifurcation event. (De)phosphorylation of MAP2 is a Michaelis-Menten function of the calcium concentration, with small changes in calcium possibly resulting in a large shift in the balance between phosphorylated and dephosphorylated MAP2. Calcium entry is assumed along the length of the dendrite (putatively through voltage-gated channels). The biochemical (de)phosphorylation pathways, involving at least CaMKII and calcineurin, are not explicitly modelled. The elongation and branching rates are given by:
where c
l
is the concentration of free tubulin at a neurite tip, B
l
is the concentration of dephosphorylated and microtubule-bound MAP2, P
l
is the concentration of phosphorylated MAP2 and P
b
, is the probability of branching.
This complex model of elongation and branching has been implemented as a coupled system of ODEs, with their numerical solution following a 'compartmental" approach, as has been used for solving the voltage equation for simulating electrical activity. In this approach, a neurite is divided into a number of short compartments, with chemical concentrations being calculated for the volume of each compartment. Chemicals move between compartments due to bulk diffusion and active transport. As a neurite elongates, new compartments are added when needed. New compartments are also created when a branching event occurs. Various strategies for when and where new compartments are added have been investigated [57]. As with PDE models, this approach also allows visualisation of concentration gradients along the length of neurites, in contrast to ODE models that only consider concentrations in the cell body and in neurite terminals [19, 20, 37, 54]. Work is needed to reconcile this compartmental approach with the numerical solution of the PDE model of neurite elongation [23, 24, 26].
