The chemotactic response of cells to graded fields of chemical cues is a complex process that requires the coordination of several intracellular activities. Fundamental steps to obtain a front vs. back differentiation in the cell are the localized distribution of internal molecules and the amplification of the external signal. The goal of this work is to develop a mathematical and computational model for the quantitative study of such phenomena in the context of axon chemotactic pathfinding in neural development. In order to perform turning decisions, axons develop front-back polarization in their distal structure, the growth cone. Starting from the recent experimental findings of the biased redistribution of receptors on the growth cone membrane, driven by the interaction with the cytoskeleton, we propose a model to investigate the significance of this process. Our main contribution is to quantitatively demonstrate that the autocatalytic loop involving receptors, cytoplasmic species and cytoskeleton is adequate to give rise to the chemotactic behavior of neural cells. We assess the fact that spatial bias in receptors is a precursory key event for chemotactic response, establishing the necessity of a tight link between upstream gradient sensing and downstream cytoskeleton dynamics. We analyze further crosslinked effects and, among others, the contribution to polarization of internal enzymatic reactions, which entail the production of molecules with a one-to-more factor. The model shows that the enzymatic efficiency of such reactions must overcome a threshold in order to give rise to a sufficient amplification, another fundamental precursory step for obtaining polarization. Eventually, we address the characteristic behavior of the attraction/repulsion of axons subjected to the same cue, providing a quantitative indicator of the parameters which more critically determine this nontrivial chemotactic response.
Growth cones are 3D hand-shaped structures, which dynamically change their conformation. Filopodia protrude from the GC membrane, continuously extending and retracting to explore the environment and to create adhesion to the substratum. These phenomena are highly complex and take place at a time frequency faster than the one we are interested into, typically seconds vs. a few minutes. For the purposes of the present model, we neglect fine local details of the shape and we represent GCs as 2D disk-like structures with diameter of. This simplification models the average shape assumed by GCs in their state prior to actual motion (the “pausing state” cited in [19]). The same assumption has been made by other authors interested in GC mathematical modelling, see [16], [33], [34]. We refer as a paradigm to the in vitro chemotactic assay for neural cells, where GCs are exposed to steady graded concentrations of a chemical ligand released by a pipette (see, e. g. , [35], [36]). Observe that this is essentially a two–dimensional situation. Polar coordinates are used, the origin of the axes being positioned in the center of the GC. The angle denotes the azimuthal coordinate on the membrane and its origin is set along the direction connecting the GC center with the pointwise source, which we always suppose to lay on the right hand side of the cell (see Fig. 1). We denote by the radius of the GC. In the model, we deal with membrane species (receptors, cyclases and calcium channels) and cytosolic species (G–proteins, cyclic AMP, cyclic GMP, kinases and their corresponding precursors). To connect membrane and cytosolic concentrations (and vice versa), we use the dimensional corrective factor (or its inverse). We use the notation to represent concentration of a certain species; moreover, we denote by the superscript the bound or active form of the molecules. The nomenclature for the species is reported in Tab. 1. In this section, we propose an extended version of the model to study the bifunctional response to a guidance cue, a phenomenon known to interest the response to netrins [27]. Engagements of receptor complexes is known to control the specificity and the polarity of the response of the neuron to the guidance cue. Here we will not deal with the details of the dynamics of the formation of the DCC–UNC5 complex, and we will always consider such a dynamics at the equilibrium. From the modelling point of view, this amounts to prescribe a priori the percentage of DCC receptors forming a complex with UNC5 (on this issue, see also the discussion in the Methods Section). Setting in the tests various percentages, we analyze a wide spectrum of situations. To describe the DCC–UNC5 complex, we follow here the idea proposed in [28], and namely that such a complex leads to cGMP synthesis, regulated by 12-hydroxyperoxyeicosatetraenoic acid (HPETE) via direct activation of guanylate cyclase [39]. Enhancement of the cGMP level causes to calcium channel closure [40]. This fact is at the origin of the significant decrease of inward calcium flux in UNC5–overexpressing GCs, repelled by netrin. Due to the substantial similarity with the sole DCC pathway, many mechanisms are modeled here in the same way. The binding reaction for the DCC–UNC5 complex (for clarity indicated just as UNC) readsThe ligated complex induces dissociation of a G–protein (probably of Rho type, see [12]), that we denote here by, givingThis process leads to the formation of HPETE. Not being yet completely explained, we consider here a second order reaction to occur (see the Methods Section for a more detailed discussion of this aspect) HPETE activates soluble guanylate–cyclase () which, in turn, catalyzes the synthesis of cGMP from the guanosine triphosphate (GTP) substrateFormation of the guanosine–dependent kinase protein (PKG) enzymatically enhances closure of calcium channels
We first demonstrate that the model correctly achieves polarization and reaches a steady state condition, where front–back differentiation is established. In Fig. 4, we plot the concentration profiles as a function of time, obtained from simulations carried out till. In Fig. 5, we plot for the different species the abscissa of the barycenter of the molecules as a function of time. Since at the initial time the distribution of molecules is homogeneous, the displacement of their spatial barycenter from represents an index of the intensity of the polarization. In Fig. 6, we plot the concentrations after as a function of. Significant polarization of the receptors takes place in tenth of minutes; polarization is inherited by all the internal species. An interesting behavior is shown by calcium channels, which undergo in all sectors a first phase of opening, reaching a fairly similar maximum value, followed by closure, more pronounced in the rear side. This mechanism might represent a sort of LEGI, global “inhibition” being constituted by collective opening and local “activation” by differential closure. Note that this is to be intended only as a qualitative interpretation (see also the discussion in [11]). The position of the barycenter of bound receptors presents a sigmoid behavior, which is characteristic of autocatalyzing processes: a first phase of relatively slow accumulation (lag time) followed by a quick growth till a steady state. This is in agreement with the experimental result of [19, Fig. 1e]. In Fig. 7, we plot the concentration profiles as a function of time for a simulation with a longer integration interval. A steady state is definitely reached. In [19], an experiment is reported, where a chelator is used to subtract the calcium available in the cytosol, obtaining a suppression in the asymmetric relocalization of receptors. To perform an in silico investigation of this experiment, we have studied the effect of the variation of the feedback coefficient. Fig. 8 (left) shows that significant receptor relocalization is obtained only above a nonzero threshold of. More importantly, the model predicts that the lack of receptor relocalization implies absence of chemotactic response and not only a weaker, but still existing, response. Fig. 8 (right) quantifies how the asymmetry in DCC receptors localization is reflected in downstream differential opening of the calcium channels. This indicates that a sufficient active relocalization of receptors is an upstream enhancing event needed to produce chemotactic response. In a second investigation, we have considered the effect of amplification due to chemical kinetics. In the signalling pathway, we may identify –loosely speaking– two families of reactions. The first family consists of stoichiometric reactions, as for example ligand–receptor or G–cyclase binding, involving one–to–one molecule synthesis. The second family consists instead of enzymatic reactions, which, involving a one–to–more molecule production, drive an internal amplification processes. We focus our attention on this latter family, taking as a representative case the cAMP production catalyzed by cyclase. We perform different simulations with a decreasing kinetic constant, which quantitatively modulates how many molecules of cAMP are produced starting from an available molecule of activated cyclase. In Fig. 9, we plot the position of the barycenter of the bound receptors as a function of. This result suggests that polarization is also indissolubly crosslinked with internal amplification. A parametric constraint appears: below a non–zero threshold in the enzymatic (that is, amplifying) efficiency of the pathway, no polarization occurs. The above investigations show that both feedback and amplification are precursor events for the achievement of polarization. Their respective actions are necessary and concurring contributions. The results discussed above were obtained setting in the model the diffusion and feedback coefficients as indicated in Tab. 2. While the diffusion coefficient can be measured as a well defined physical parameter, it is much more difficult to quantify the feedback parameter. This fact has a strong implication, since the ratio of the two parameters influences the dominating behavior of the system. To fix the ideas, we perform the tests varying the feedback coefficient, and we keep constant the diffusion coefficients, (that we denote here for simplicity as). The results of the simulations are shown in Fig. 10, where the abscissa of the barycenter (in) of the bound receptors is plotted as a function of. They suggest that under a certain threshold, diffusion overwhelms drift, leading to an unbiased receptor distribution on the membrane, that is, as if the external field were uniform. We use the model to study the time scales that characterize the process on the front and back sides, respectively. To perform this mathematical analysis, we consider a simplified version of the DCC model, neglecting diffusion and feedback terms. By doing so, we yield a system of ordinary differential equations, decoupled sector by sector. In this study, we prescribe a-priori an asymmetric receptor distribution to describe the polarized situation reached after a sufficient time of exposure to the cue. In particular, we start from the steady state distribution of receptors obtained from the simulation of the DCC model with. We compute in each sector the eigenvalues of the Jacobian matrix of the system in correspondence to its steady state. All the eigenvalues are real negative, indicating that the steady state is an attractive point. Based on the principal component of the corresponding eigenvector, we associate a chemical species with each eigenvalue. Then, using standard tracking techniques [41], we follow the variation of each eigenvalue along the GC perimeter. In Fig. 11, we plot the modulus of the eigenvalue associated with each species as a function of the angle. The eigenvalue associated with the slowest process on the front side () appears to be connected to PKA, while on the back side () it appears to be connected to. Observe that all the eigenvalues undergo a variation along the angle, even if for some of them this is not apparent in the logarithmic scale, required to appreciate the different relative behaviors. The graph shows the strong variation of the eigenvalue connected with. Moreover, the general trend of reduction of the absolute values passing from the front side to the back side indicates that the front dynamics is faster than the back one. We use the coupled DCC–UNC5 model to study the response of the system to the (bifunctional) netrin cue when the DCC–UNC5 complex is formed. Denoting by #DCC and #UNC the number of receptors on the membrane belonging to the two populations, respectively, we introduce the quantitywhich represents the fraction of DCC–UNC5 receptors. Neurons which display a chemoattractive response when exposed to a netrin gradient are typically characterized by low values of. Neurons genetically manipulated to overexpress the UNC5–type receptors, which display a chemorepulsive response when exposed to the same netrin gradient, are typically characterized by values of near. As a representative situation of the first case, we set. In Fig. 12, we plot the concentration contours of the main species as a function of angular position and time. This model suggests the following explanation of the attractive behavior: DCC receptors migrate toward the source, while DCC–UNC5 ones migrate away from it, causing a differential opening of calcium channels on front vs. back. Then, we consider the dual case, setting. In Fig. 13, we plot the concentration contours of the main species as a function of angular position and time. The model suggests a dual behavior with respect to the situation with: DCC–UNC5 receptors migrate toward the front, while DCC receptors migrate toward the back, giving rise to the repulsive response.
We have proposed a mathematical model to study the polarization phenomenon triggered by the exposure of GCs to an external cue, taking as a paradigm the in vitro chemotactic assay. The key hypothesis is that symmetry breaking occurs as early as at the level of transmembrane receptors, which undergo a biased distribution after exposure to the cue. This finding appears in the recent work [19], (see also the analysis of [42]), where experimental evidence of such a process is provided and an idea is proposed of an autocatalyzing loop connecting receptors and downstream actin dynamics. Our main contribution is to quantitatively demonstrate via the mathematical model how such a loop is able to achieve GC “front–back” polarization. More precisely, we assess the fact that spatial bias in receptors is a precursor necessary event for chemotactic response, so that upstream gradient sensing and downstream cytoskeleton dynamics cannot be decoupled. Moreover, we analyze further crosslinked effects and, namely, the contribution to polarization of internal enzymatic reactions, which entail activation of a one–to–more production of molecules. The model shows that the enzymatic efficiency of such reactions must overcome a threshold in order to produce a sufficient amplification, which is another fundamental precursor step for obtaining polarization. A simplified version of the model is used to provide preliminary indications about time scales in the front and back processes, via eigenvalue tracking. The findings suggest that the two processes take place with different speeds, the front one being faster than the rear one. This element could play a significant role, even if further investigations should be carried out, A preliminary analysis of the role of diffusion vs. convection is also sketched, establishing the nature of a strongly advection–dominated system. Eventually, we have proposed an extension of the model to address a peculiar behavior arising when two families of receptors interact to produce the response to the same cue. The study of this case allows to propose some ideas on the mechanism of chemorepulsion, as the synergistic interaction of pathways, that contribute to better understand this less studied phenomenon. The present model demands for improvement. We have run simulations considering a time dependent gradient of ligand. The laboratory experiments (among the very few on this issue for neural cells) reported in [19, Fig. 1g], show reversibility of the distribution of the receptors, upon removal of the gradient. Our model tested under these conditions tends instead to maintain the receptor asymmetric configuration and removal of the cue itself is demanded. This fact has led us to conjecture the need for introducing inhibitory effects on receptors. At present, we have indeed assumed that the GC has a pool of receptors capable at each time to bind the ligand and we have not considered their physiological cycle, comprising deactivation by phosphorylation and recycling back in the membrane (as done for example in [29]). Attempts have been made to consider such processes, establishing a life–time for each receptor, after which the receptor is removed and substituted by another one without bias in its angular position on the membrane (keeping in this way the total number of receptors constant in time). Including these effects in the model does not seem to significantly change the results, but deeper investigations on these ideas should be carried out, possibly supported by the availability of further information from the biochemical viewpoint (as done in the context of nonneural cells, for example, in [43]).
To derive initial conditions for the mathematical model, which are not readily available from experiments, we start assuming the total number of molecules per cell (indicated with the subscript) of the following species to be conserved during the integration time, leading to, for all (7a) (7b) (7c) (7d) (7e) (7f) (7g) Total values used in the simulations are reported in Tab. 3. Moreover, we assume that: Using assumptions (1), (2) and (3), we recover from the total values the corresponding concentrations at the left hand sides of Eqs. (7a–g) at. Values are in Tab. 4. As for the boundary conditions, since the GC is modeled as a circular structure, concentrations and fluxes must coincide at and. Thus, for each species, we impose the periodic boundary conditions Several kinetic constants entering the model are available from literature references for the same reactions we are dealing with or for very similar reactions. The value of the other parameters has been estimated, based on the following considerations: Tab. 5 summarizes the value of the kinetic constants used in the model. We conclude this discussion exploring the overall influence of the parameters over the model predictions. To do this, we consider the DCC model and we perform a set of 300 trial tests prescribing a random variation of all kinetic constants. We monitor the displacement of the barycenter of the bound DCC receptors after 2 h. The results are shown in Fig. 15, where on the axis we report the total parameter variation computed as, being the total number of parameters, the value of the parameter and its perturbed value. A percentage of the tests show a perturbation in the barycenter displacement lower than the. Both DCC and DCC–UNC5 models constitute nonlinear time dependent diffusion–advection–reaction systems of partial differential equations. Their numerical approximation is a very challenging task: to start with, we partition the GC perimeter into angular sectors (see Fig. 1) and we discretize the spatial derivative operators using finite differences, with a node collocated at the center of each sector. Observe that conservation relations like Eqs. (7) do not hold sector by sector, but they rather apply to the integral on the angle (which represents the number of molecules per cell). Attention must be paid to the fact that the feedback terms in Eqs. (4a) and (5a) dominate the diffusive terms. To avoid spurious oscillations, upwind finite differences or a sufficiently fine discretization should be adopted [45]. Once finite difference discretization is carried out, a system of coupled nonlinear ordinary differential equations is obtained. Due to the different speeds in the reaction dynamics (refer for this issue also to the Results and Discussion), these systems are very stiff and require an implicit time integrator. We have adopted the ode15s Matlab routine with adaptive choice of the time integration step. Numerical evaluation of the Jacobian matrix has been used for the linearization. The Matlab software package developed by the authors can be made available upon request.
The ability of cells to respond to chemical signals present in the environment is of upmost importance for life. In the developing embryo, cells crawl along graded fields of chemical cues to aggregate into organized patterns. This process is an example of chemotaxis. It is a complex phenomenon, where external signals are transduced into internal chemical pathways leading to directional movement. Differential reorganization of the internal structure is called polarization, and it involves regulatory proteins as well as cytoskeletal elements. In this work, we propose a mathematical and computational model for the quantitative study of chemotactic pathfinding in neural cells. Our starting point is the recent finding that, for such cells, an early polarization event is the redistribution on the membrane of cue–ligated receptors, transported by the cytoskeletal structures, which act as a sort of conveyor belt. We show that this proposed mechanism, connecting in a closed loop cue sensing and cytoskeleton dynamics, is qualitatively and quantitatively adequate to produce polarization. We also investigate the role of the internal biochemical chain in producing signal amplification and its tight interlacing with polarization. An extension of the model is used to study chemotactic behaviors as the attractive/repulsive response of axons exposed to the same cue.
computer science/applications
cell biology/cell signaling
mathematics
neuroscience/neurodevelopment
computational biology
chemical biology/chemical biology of the cell
Autocatalytic Loop, Amplification and Diffusion: A Mathematical and Computational Model of Cell Polarization in Neural Chemotaxis