Modeling material transport regulation and traffic jam in neurons using PDE-constrained optimization.

The intracellular transport process plays an important role in delivering essential materials throughout branched geometries of neurons for their survival and function. Many neurodegenerative diseases have been associated with the disruption of transport. Therefore, it is essential to study how neurons control the transport process to localize materials to necessary locations. Here, we develop a novel optimization model to simulate the traffic regulation mechanism of material transport in complex geometries of neurons. The transport is controlled to avoid traffic jam of materials by minimizing a pre-defined objective function. The optimization subjects to a set of partial differential equation (PDE) constraints that describe the material transport process based on a macroscopic molecular-motor-assisted transport model of intracellular particles. The proposed PDE-constrained optimization model is solved in complex tree structures by using isogeometric analysis (IGA). Different simulation parameters are used to introduce traffic jams and study how neurons handle the transport issue. Specifically, we successfully model and explain the traffic jam caused by reduced number of microtubules (MTs) and MT swirls. In summary, our model effectively simulates the material transport process in healthy neurons and also explains the formation of a traffic jam in abnormal neurons. Our results demonstrate that both geometry and MT structure play important roles in achieving an optimal transport process in neuron.

Introduction

The neuron exhibits a highly polarized structure that typically consists of a single long axon and multiple dendrites which are both extended from its cell body. Since most of the materials necessary for the neuron are synthesized in the cell body, they need to experience long-distance transport in axons or dendrites to reach their effective location[1,2]. The intracellular material transport is therefore especially crucial to ensure necessary materials are delivered to the right locations for the development, function, and survival of neuron cells. The disruption of intracellular transport can lead to the abnormal accumulations of certain cellular material and extreme swelling of the axon, which have been observed in many neurological and neurodegenerative diseases such as Huntington’s, Parkinson’s, and Alzheimer’s disease[3-7]. Therefore, it is essential to study and understand mechanisms of the transport function and dysfunction.

Recent studies have shown that the neuron is critically dependent on molecular motors to transport various materials along the longitudinal cytoskeletal structure like microtubules (MTs)[8-10]. MTs are long and polarized polymers with biophysically distinct plus and minus ends[11-13]. The polarity of MTs can decide the preferred direction in which an individual molecular motor moves. For instance, molecular motors from the kinesin and dynein superfamilies have been identified to convey materials along MTs towards their plus and minus ends respectively[14,15]. Inspired by these findings, there have been many mathematical models proposed to quantitatively study the motor-driven transport process and understand the pathology of neuron diseases. For instance, the partial differential equations (PDEs) of linear reaction-hyperbolic form have been used to approximate the traveling waves of a single moving species[16]. This model was further extended to account for multiple moving species[17] and their diffusion[18,19]. Based on PDE-based transport, stochastic models have also been developed for both axonal transport[20,21] and dendritic transport[22,23]. In addition, several mathematical models were developed to simulate material transport in unhealthy neurons. Xue et al. presented a stochastic model to explain the segregation of MTs and neurofilaments in neurological diseases[24]. Bertsch et al. proposed to couple Smoluchowski equations and kinetic-type transport equations to study the onset and progression of Alzheimer’s disease[25].

Though the aforementioned PDE and stochastic models can successfully simulate and explain certain phenomena during transport, most of these models were solved only in simple one-dimensional (1D) or 2D domains without considering the complex neuron geometry. Recent developments in numerical methods allow us to obtain accurate solution of PDEs in complex geometries. Specifically, isogeometric analysis (IGA)[26] directly integrates geometric modeling with numerical simulation and achieves better accuracy and robustness compared to the conventional finite element method (FEM), making it a perfect tool to tackle the highly branched neuron geometry. In particular, IGA performs simulation with different types of splines as basis functions instead of Lagrange polynomials used in conventional FEM. The same smooth spline basis functions[27] used for both geometrical modeling and numerical simulation lead to accurate geometry representation with high-order continuity and superior numerical accuracy in simulation. Therefore, IGA has been extensively used in shell analysis[28-31], cardiovascular modeling[32-37], neuroscience simulation[38,39], fluid-structure interaction[40-43], as well as industrial application[44,45]. Truncated T-splines[46,47] were developed to support local refinement over unstructured quadrilateral and hexahedral meshes. Blended B-splines[48] and Catmull-Clark subdivision basis functions[49] were investigated to enable improved or even optimal convergence rates for IGA. With the advances in IGA, we developed an IGA-based simulation platform to accurately reconstruct complex neuron geometries and solved a 3D motor-assisted transport model within them[38]. We also developed a deep learning framework based on the IGA simulation platform to predict the material transport process in complex neurite networks[50]. The results from our IGA solver showed how the complex neuron geometry affects the spatiotemporal material distribution at neurite junctions and within different branches. However, the motor-assisted model only provides a simplified model of the actual transport process but ignores the active regulation from neuron itself.

To model the active regulation from neurons to control the transport process, we propose to use PDE-constrained optimization (PDE-CO). PDEs are commonly used in science and engineering to mathematically represent biological and physical phenomena. Recent advances in numerical methods and high-performance computing equip the development of large-scale PDE solvers. As a result, PDE-CO problems arise in a variety of applications including optimal design[51-53], optimal control[54-56], and inverse problem[57,58]. In particular, PDE-CO has important biomedical applications in exploiting valuable information from real medical data. For instance, Hogea et al. presented a PDE-CO framework for modeling gliomas growth and their mass-effect on the surrounding brain tissue[59]. Kim et al. proposed a transport-theory-based PDE-constrained multispectral imaging algorithm to reconstruct the spatial distribution of chromophores in tissue[60]. Melani utilized the blood flow data and solved a PDE-CO problem based on fluid-structure interaction to estimate the compliance of arterial walls in vascular networks[61]. PDE-CO problems was also used to model tumor growth model by fitting the numerical solution with real experiment data and estimating unknown parameters in the model[62,63].

In this study, we develop a novel IGA-based PDE-CO framework that effectively simulates the material transport regulation and investigates the formation of traffic jams and swirl during the transport process in complex neurite structures. Our simulation reveals that the molecular motors and MT structure play fundamental roles in controlling the delivery of material by mediating the transport velocity on MTs. The defective transport on MTs can cause material accumulation in a local region which may further lead to the degeneration of neuron cells. Combined with geometry of the neurite network, the motor-assisted transport on MTs controls the routing of material transport at junctions of neurite branches and distributes transported materials throughout the networks. Therefore, our study provides key insights into how material transport in neurite networks is mediated by MTs and their complex geometry. In summary, there are three main contributions in this paper.

Our PDE-CO model introduces a new objective function to simulate two transport control mechanisms for (1) mediating the transport velocity field; and (2) avoiding the traffic jam caused by local material accumulation. The control strength can be adjusted through two penalty parameters in the objective function and the impact of these parameters is also studied;

Our model introduces new simulation parameters to describe the spatial distribution of MTs, which enables the simulation of traffic jam caused by abnormal MT structure such as MT reduction and MT swirls; and

Our study develops an IGA optimization framework for solving the PDE-CO problem in complex neuron geometries. The optimization framework is transformative and can be extended to solve other PDE-CO models of cellular processes in complex neurite networks.

IGA-based material transport optimization in neurons

Our interest lies in the transport of particles along an axon or dendrite in neuron cell. In our previous work, we simulated the material transport process using a macroscopic molecular-motor-assisted transport model without any transport control[38]. Built upon this transport model, we propose a novel transport optimization model to further study the transport control mechanism of neuron and predict the formation of a traffic jam in abnormal neurons. Assuming the open set ( or 3) represents the d-dimensional internal space of the neuron, , and are referred as the “state variables” while are referred as the “control variables”. The proposed optimization problem is described as where is a predefined velocity field inside neuron; , and are the spatial concentrations of free, incoming (relative to the cell body; retrograde), and outgoing (anterograde) particles, respectively; is the diffusion coefficient of incoming and outgoing materials; and are velocities of incoming and outgoing particles, respectively; and are rates of MT attachment and detachment of incoming and outgoing materials, respectively; represents the density of MTs used for motor-assisted transport; represents the control forces that mediate the material transport; is viscosity of traffic flow; , represent the degree of loading at inlet and outlet ends, respectively[18]; and , represent the boundary value of at inlet and outlet ends, respectively. On the right hand side of Eq. (1d), the concentration gradient term accounts for the porosity induced by the material, which models the behavior of molecular motors on MT to speed up when the material concentration decreases while slow down when increases. In this study, we assume the MT system is unipolar that leads to a unidirectional material transport process and ignore , , , , terms in Eq. (1b)–(1f). The selection of parameters is based on the study of Smith and Simmons[18]: the Einstein-Stokes relation for a 1 m sphere in water gives m/s and we reduce the value to m/s considering the irregular neuron topology and a bigger cytoplasmic viscosity; the attachment rate is set to be and the detachment rate is set to be . The default values of simulation parameters are summarized in Table 1.

Table 1

Simulation parameters utilized in computations.

Parameter	Description [Unit]	Default value
D_±	Diffusion coefficient of incoming and outgoing materials [μm^2/s]	0.1
k_±	Attachment rate to the MTs that transport materials in the positive (+) and negative (−) directions [s^-1]	1.0
k±′	Detachment rate from MTs for materials that move in the positive (+) and negative (−) directions [/s]	0.1
l_±	Density of MTs used for motor-assisted transport	1.0
μ	Viscosity of the traffic flow [μ/m/s]	0.1
λ_i	Degree of loading at inlet end	2.0
λ_o	Degree of loading at outlet end	2.0
n_i	Boundary value of n_0 at inlet end [mol/μ/m2]	1.0
n_o	Boundary value of n_0 at outlet end [mol/μ/m2]	0.0
α	Penalty parameter for the cost to control high concentration gradient	1.0
β	Penalty parameter for the cost of control force f_±	1.0

Herein, we account for active regulation from neuron in the objective function (Eq. 1a), and we assume the optimal material transportation within neuron can be achieved by solving the proposed optimization model. The first term in Eq. (1a) measures the difference between and the predefined optimal velocity field . It serves as a velocity control mechanism that neuron expects to achieve the predefined velocity field during transport. The second term measures the cost from concentration gradient within the entire neuron cell. It serves as a traffic jam control mechanism that the neuron can improve local traffic jam by detecting and avoiding high concentration gradient in the entire geometry. The value of parameter represents to what extent we want to optimize the transport process and avoid traffic jams. The third term is a regularization that measures the control forces applied by neuron to mediate the transport. The value of parameter represents how much the neuron can affect the transport velocity. To introduce traffic jams in neurons, we modify the simulation parameters in the governing equations. In this study, we modify the spatial distribution of to model the traffic jam caused by abnormal MTs such as the reduction of MTs and MT swirls during transport.

We employ the “all-at-once” method[64,65] and IGA to formulate and solve the optimization model (Eq. 1a) with PDE constraints (Eqs. 1b–1f) simultaneously. We first discretize the objective function to obtain its approximation . We also discretize PDE constraints (Eqs. 1b–1f) to obtain their weak form . Then, we build a discrete Lagrangianwhere is the Lagrange multiplier and is also referred to as the “adjoint variable”. By taking derivatives of the discrete Lagrangian with respect to state, control, and adjoint variables and setting the resulting expressions to zero, we obtain the first-order conditions, or Karush-Kuhn-Tucker (KKT) conditions. The resulting KKT system is then solved using the GMRES[66] solver implemented in PETSc[67]. In this study, we focus on solving the proposed optimization model in 2D neuron geometries. To handle the ill-conditioning issue of the KKT system, we implemented the preconditioner following the all-at-once method[64]. The stopping tolerance for the optimization problem is set to be . We do not find other local minima when solving the optimization problem. However, due to the nonlinear constraints in our model, it is possible to have other local minima. Overall, there are 7 unknowns in the PDE-CO model, including , , , , , and . In particular, , and are ignored when no MT swirl is introduced. The computation usually needs 1 to 3 CPU Nodes (128 cores per node) and takes 8-12 hours to finish.

As shown in Fig. 1, we use a bifurcation example to illustrate the pipeline of our simulation. We first generate a control mesh and reconstruct the neuron geometry with Truncated Hierarchical B-splines (THB-spline) by utilizing the geometry information stored in a SWC file. The SWC file is widely used to store neuron morphologies including vertices and the associated diameters on the skeleton of the neuron. We can obtain the SWC files for various real neuron geometries from the NeuroMorpho database[68]. The raw SWC file needs to be pre-processed to ensure no duplicated vertices or overlapping skeleton exist in the geometry. During the geometric modeling of our workflow, we take the cleaned-up neuron skeleton as input and use the skeleton-based sweeping method[32] to generate quadrilateral control mesh of the neuron geometry. Then, we build THB-spline on the quadrilateral mesh[30,31] for the final representation of the neuron geometry. Once the spline information for the geometry is obtained, we run a steady-state Navier–Stokes solver to generate the pre-defined velocity for the optimization. We then use the default simulation parameters in Table 1 and modify the spatial distribution of in the red circle regions to introduce traffic jam. Finally, we run the optimization solver and obtain the velocity field and concentration distribution. In this paper, we apply the pipeline to various neural structures with material transport regulation, traffic jam and MT swirl. All simulations are conducted on the XSEDE (Extreme Science and Engineering Discovery Environment) supercomputer Bridges at the Pittsburgh Supercomputer Center[69,70].

Figure 1

An overview of the material transport control simulation in a bifurcation geometry. The traffic jam is introduced by reducing MTs in the red dashed circle region. Color bars unit for velocity field: m/s and concentration: mol/m.

Results

Simulation of material transport regulation and traffic jam

We first simulate the normal material transport and the abnormal transport with traffic jam in a single pipe geometry (Fig. 2). The predefined velocity field for both cases is computed by solving a steady-state Navier–Stokes equation and the result is shown in Fig. 2A. The other simulation parameter settings are summarized in Table 1. The computed velocity field and the distribution of concentration in the normal transport are shown in Fig. 2B,E. To model traffic jam caused by the reduction of MTs, the distribution of along the pipe is defined as shown in Fig. 2D. The velocity field and material distribution results in the abnormal transport are shown in Fig. 2C,F. The convergence curve for the gradient of objective functions is shown in Supplementary Fig. S6 to verify the result. The comparison between normal and abnormal transport shows that the velocity magnitude decreases in the red dashed circle region due to the reduced number of MTs, and this further leads to accumulation of the material in this area.

Figure 2

Simulation of material transport and parameter analysis in a single pipe geometry. (A) The predefined velocity field . Black arrow points to the inlet of the pipe. The computed velocity field in (B) a healthy neuron and (C) an abnormal neuron with reduced MTs in the red dashed circle region. (D) Distribution of to model the traffic jam caused by the reduction of MTs. Distribution of concentration () in (E) a healthy neuron and (F) an abnormal neuron with reduced MTs in the red dashed circle region. (G–I) The concentration curve () on the centerline of the single pipe affected by different settings of (G) ; (H) ; and (I) . Unit for color bars: (A–C) m/s and (E, F) mol/m.

As shown in Fig. 2G–I and Fig. S1, we also perform parameter analysis using the single pipe geometry to study the influence of simulation parameters on the material distribution results. In particular, we focus on three parameters that may have significant effect when dealing with traffic jam caused by the reduction of MTs. The values selected for these parameters are displayed in Table 2. We assume the active regulation from neuron is less dominant than natural transport via diffusion or MTs, and thus select two smaller values for and compared to the default values in Table 1. Regarding the value selection of , we refer to the values utilized in[71] and ensure the selected values stay within a biologically realistic range. Figure 2G shows the effect of the penalty parameter of the concentration gradient cost, , on the concentration distribution. One can see that the decrease of leads to a severer material accumulation around the region with reduced MTs in the single pipe geometry. We also find that the concentration gradient becomes larger around the traffic jam region, which indicates that there is less control over the concentration gradient due to the decrease of .

Table 2

Value selection for parameter study.

Parameter	Value selection
α	1, 0.1, 0.01
β	1, 0.1, 0.01
k_±/k±′	Fix k_±=1.0s^-1, let k_±/k±′ = 1, 10, 100

Figure 2H is similar to Fig. 2G but shows the effect of the penalty parameter of the control force, , on the concentration distribution. We find similar phenomena that the traffic jam gets worse when decreases. By comparing Fig. 2G with 2H and Fig. S1A with S1B, we find has a greater influence on the concentration than when decreasing both parameters by the same amount. Since affects the control force in Eq. (1d) while affects the concentration in Eqs. (1b) and (1c), the result indicates that the regulation of transport velocity on MTs is vital to achieve the optimal material transport process in neuron. In addition, we study the impact of and on the distribution of the control force f and the resulting velocity profile as shown in Supplementary Fig. S2. When increases, the velocity increases due to the decrease of . f also decreases since less control is needed to regulate the transport. When increases, f becomes smaller and the overall velocity decreases. We also notice that with larger , the distribution of f is stabler and the sharp velocity change is mitigated in the traffic jam region, which may explain the reduction of . Though the increase of regularization reduces the overall transport, the material accumulation is alleviated and has less hazard to the neuron, which is an expected control from neuron.

Figure 2I shows the effect of the ratio between the attachment rate and detachment rate, , on the material concentration. We find that when increases, the location of maximum concentration moves toward right, which indicates the decrease of detachment rate causes more material get attached to MTs and transport faster as expected in[18]. However, the reduction of MTs slows down the motor-assisted transport on MTs and results in worse traffic jam. Interestingly, when decreases from 10 to 1, we also observe a similar traffic jam phenomena. The possible reason is that the increase of causes more material transported via free diffusion. Although free diffusion helps to transport the material farther along the branch, the slow diffusion speed limits its ability to mitigate the traffic jam caused by the reduction of MTs.

To account for morphological effect on the transport process, we simulate the normal material transport and the abnormal transport with traffic jam in two neuron tree structures as shown in Figs. 3 and 4. The predefined velocity fields for both geometries are shown in Figs. 3A and 4A. To quantitatively study the influence of traffic jam on the material concentration among tree structures, we also plot the concentration distribution curves along the centerline from the inlet to each outlet of these two neurons. In each curve plot, we compare the distribution between the normal transport and the abnormal transport with traffic jam, as shown in Figs. 3E and 4E. For both cases, we model traffic jam by reducing the number of MTs () used for transport in the red dashed circle regions. As a result, a sudden decrease of velocity (Figs. 3C and 4C) and material accumulation (Figs. 3E and 4E) can be observed in these regions. By observing the distribution curve of the outlets downstream the traffic jam region (curve plots 1–4 of Fig. 3E and 3–8 of Fig. 4E we find that the reduced number of MTs not only causes high concentration in the local region, but also decreases the material concentration along the downstream of traffic jam region. The distribution curves of the other outlets (curve plots 5 of Fig. 3E and 1, 2, 9, 10 of Fig. 4E) demonstrate that more materials are transported to these outlets to minimize the hazard of traffic jam. The result also shows that materials rely on motor-assisted transport in longer branches of neurons and the directional transport on MTs contributes significantly to the entire transport process.

Figure 3

Simulation of material transport and parameter analysis in a neuron tree extracted from NMO_54504. (A) The predefined velocity field . Black arrow points to the inlet of the neuron tree. The computed velocity field in (B) a healthy neuron and (C) an abnormal neuron with reduced MTs in the red dashed circle region. Distribution of concentration () in (D) a healthy neuron and (E) an abnormal neuron. We also compare the concentration curve on the centerline from the inlet to every outlet between normal and abnormal transport in (E). The red dashed curve shows the centerline from the inlet to one of the outlets and each outlet is indexed by a unique number. (F–H) The concentration curve () on the centerline from inlet to outlet 2 affected by different settings of (F) ; (G) ; and (H) . Unit for color bars: (A–C) m/s and (D, E) mol/m.

Figure 4

Simulation of material transport in a neuron tree extracted from NMO_54499. (A) The predefined velocity field . Black arrow points to the inlet of the material. The computed velocity field in (B) a healthy neuron and (C) an abnormal neuron with reduced MTs in the red dashed circle region. Distribution of concentration () and the concentration curve on the centerline of the circled region in (D) a healthy neuron and (E) an abnormal neuron. We also compare the concentration curve on the centerline from the inlet to every outlet between normal and abnormal transport in (E). The red dashed curve shows the centerline from the inlet to one of the outlets and each outlet is indexed by a unique number. Unit for color bars: (A–C) m/s and (D,E) mol/m.

As shown in Fig. 3F,H and Figs. S3 and S4, we also perform parameter analysis on the concentration distribution in the neuron tree structure. Similar to the parameter analysis in single pipe geometry, we study the influence of three parameters on the concentration distribution and the selected values are listed in Table 2. To quantitatively study the influence, we also plot and compare the concentration curves on the centerline from inlet to outlet 2 of the neuron tree. We obtain similar results as in single pipe geometry that the decrease of or leads to a severer material accumulation around the region with reduced MTs, and shows greater effect than on the concentration distribution. In addition, we observe in Figs. 3E,F and S3 that when or increases, the concentration in the bottom long branch (pipe 5) has slight increase, indicating that more material is transported to the unaffected region to mitigate the traffic jam in other branches. As shown in Supplementary Fig. S4, the increase of regularization parameter reduces the control force f. The overall velocity is also reduced but we observe that velocity in pipe 5 keeps increasing and is close to the predefined value when , which explains the slight increase of the concentration in pipe 5. In Fig. 3H, we also find that the maximum concentration location moves downstream slightly when increases, and either increasing or decreasing intensifies the traffic jam.

Simulation of traffic jam with MT swirls and local swelling

Recent studies have shown that the formation of MT swirls can lead to accumulation of transported material and cause local swelling of neuron geometries[72]. In our model, we modify the spatial distribution of and enlarge the radius of neuron in a local region to simulate the effect of MT swirls and local swelling on the transport process. We explain the simulation setting by using a straight pipe geometry with MT swirls and swelling in the middle region, as shown in Fig. 5. We assume the normal transport is unidirectional from left to right ( direction, red arrow in Fig. 5A). Due to the MT swirls in the middle region, the transport path is extended by two segments: one segment reverses to transport the material from right to left (− direction, blue arrow in Fig. 5A) and the other segment transports in the normal direction from left to right. Therefore, we increase the values of and along the longitudinal direction in the swelling region to describe the transport path change caused by swirling. We also assume that the swirl direction is counter-clockwise and assign an asymmetric distribution of on the cross-section in the swelling region. In particular, is higher on the bottom of cross-section while is higher on the top of cross-section, as shown in Fig. 5A. We perform simulation with the new parameter setting and compare with the results of normal transport in the same geometry. The velocity field and concentration distribution of normal and abnormal transport are compared in Fig. 5B,C, respectively. The decrease of velocity and material accumulation can be observed in the swollen region. We also find that the velocity magnitude is not symmetric anymore due to the MT swirls in abnormal transport. In Fig. 5D, we plot the velocity streamline with concentration distribution in the zoomed-in swollen region for both normal and abnormal transport. Compared to the uniform velocity streamline in normal transport, the velocity displays vortex pattern in the abnormal transport, which reflects a longer transport distance due to MT swirls. We also find that the swirl of velocity streamline usually happens in the high concentration region, which implies that the material accumulation is caused by the vortex-shape velocity field. Moreover, we study the effect of regularization parameters on transport and the results are shown in Fig. 5 and Supplementary Fig. S5. In Fig. 5E–J, we compare the concentration distribution obtained with different and along three different curves labeled in Fig. 5C. We observe that the concentration is the highest in Curve 3 while the lowest in Curve 1, which is caused by non-uniform distribution of and on the cross section. We also find that most curves have two peaks that occur at the two ends of the swelling region. At the left end of the swelling region, the sudden increase of makes some materials to move left along “−” direction MT and thus causes traffic jam. At the right end of the swelling region, the decrease of both and causes more materials to get detached from MT and the slow diffusion causes the second material accumulation. By comparing Fig. 5E–J, we find that the increase of reduces the maximum concentration value and the concentration gradient, while the increase of reduces the range of high concentration region. As shown in Fig. S5, the overall transport velocity increases when increases and the transport velocity distribution becomes stabler when increases, which can explain the improvement of traffic jam.

Figure 5

Simulation of material transport in a straight pipe with swelling in the middle region. (A) The simulation setting for modeling MT swirls. The red and blue arrows show the transport path along swirly MTs. Due to the MT swirls in the region, both and are increased along centerline and their distributions on cross-section are also modified. (B,C) The computed velocity field and concentration distribution in the swollen geometry. Three different curves are labeled for concentration plot. (D) The velocity streamline and concentration distribution in the swollen region. Different color maps are used to distinguish between velocity and concentration. (E–J) The concentration curve () along three different curves from inlet to outlet affected by different settings of (E–G) and (H–J) . Unit for color bars: Concentration: mol/m; Velocity: m/s.

As shown in Fig. 6, we then apply the same approach to simulate the normal and abnormal transport with MT swirls in two neuron tree structures with local swelling. The swelling is introduced by increasing the skeleton radius in the red dashed circle regions. We also assume a counter-clockwise MT swirl in these swollen regions and modify the distribution of accordingly. For each model, we simulate the abnormal transport process due to MT swirls to obtain velocity field and concentration distribution results and compare with the result of normal transport in the same geometry. By comparing Fig. 6A,C with Fig. 6B,D, we find the velocity magnitude decreases and material accumulates in the swollen region. In other branches that are not downstream the swollen region, the material concentration also increases to mitigate the traffic jam in the swollen region. In addition, similar to the results in straight pipe with swelling geometry (Fig. 5D), we also observe that the velocity streamline with vortex pattern matches with the high concentration region (Fig. 6B,D). These results illustrate that the MT swirls lead to the circular transport velocity field in a local region which not only extends the transport distance but also traps the material and causes traffic jam.

Figure 6

Simulation of material transport in neuron trees extracted from (A,B) NMO_54504 and (C,D) NMO_54499 with swelling in the red circle regions. The first column shows the computed velocity field and black arrow points to the inlet of the material. The second column shows the concentration distribution. The last column shows the velocity streamline and concentration distribution in the swollen region. Different color maps are used to distinguish between velocity and concentration. Unit for color bars: Concentration: mol/m; Velocity: m/s.

Discussion

In this paper, we develop a PDE-constrained optimization model to simulate material transport control in neurons. Using our simulation, we examine both normal and abnormal transport processes in different geometries and discover several spatial patterns of the transport process. Our results show the formation of traffic jams due to the reduction of MTs and MT swirls in the local region. We also observe how the traffic jam affects the spatial patterns of transport velocities that in turn drives the transported materials distributed distinctly in different regions of neurite networks to mitigate traffic jam. By solving the proposed new optimization problem, we build a more realistic transport model for neurons by including active traffic regulation. In particular, we assume the active regulation to be an optimization process and include the potential active regulation mechanism into the objective function of the PDE-CO model. Though it is challenging and time consuming to solve the inverse problem, the model can determine a more plausible distribution of the transport control variables within neuron and provide a more reliable explanation for the active regulation mechanism. Moreover, since the objective function measures the difference between the desired and simulated results, we can further integrate the experimental or clinical data with this model to provide more realistic simulations. For instance, we can use the velocity from clinical data as input to approximate the control variables and find the potential abnormal region in the neuron that causes the disease. Herein, the model is successfully applied to complex 2D neuron geometries and provides key insights into how neuron mediates the material transport inside its complex geometry.

Our study shows that MTs have a major impact on the material transport velocity and further affect the material concentration distribution. As shown in Fig. 2, the reduction of MTs in the middle of the single pipe slows down the transport velocity downstream and leads to traffic jam in the middle region. When the neuron has more branches in its geometry (Figs. 3 and 4), the reduction of MTs in one branch has a similar influence on the transport downstream the branch. However, we observe an increase in transport velocity and material concentration in other branches, indicating that the active regulation from neuron takes effect to avoid traffic jams. In addition, we perform parameter analysis to study the influence of different simulation parameters on the material concentration distribution. The ratio between the attachment rate k and detachment rate affects the amount of material transported via MTs or free diffusion. This will affect the overall transport speed and material distribution due to the different transport behavior between motor-assisted transport and free diffusion. The penalty parameters and affect the ability of neuron to handle traffic jams. has a greater influence on the traffic regulation compared to since it directly affects the transport velocity on MTs, this again verifies the vital role of MTs during the intracellular transport process. Our model can also model the influence of diverse neuron topologies on material distribution. For the transport in healthy neurons (Figs. 3,4B,D), the magnitude of transport velocity is different among branches due to the asymmetric geometry. The different velocity magnitude further contributes to the distinct material concentration in different branches. In particular, we find that shorter branches tend to have faster transport speed and higher material concentration, which may result from the high demand of materials for their growth.

Our study also successfully simulates and provides reasonable explanation on the traffic jam caused by MT swirls. We assume the counter-clockwise MT swirls exist in a local region of neuron geometry which cause traffic jam and geometry swelling. The spatial distribution of MT density () and neuron geometry are modified accordingly to model this phenomena. We compare the simulation result of abnormal transport on swirly MTs with normal transport and find that MT swirls have severe impact on the transport velocity field. Compared to the uniform velocity streamline in normal transport, the abnormal transport exhibits a streamline with counter-clockwise vortex pattern (Figs. 5D, 6B,D), which is caused by the counter-clockwise MT swirls. This circular streamline not only extends the transport distance but also traps the material in the local region, and therefore explains why high concentration region matches with the circular streamline pattern.

Our study develops an IGA solver (available at https://github.com/CMU-CBML/NeuronTransportOptimization) for solving the PDE-CO problem in complex neuron geometries. Specifically, we adopt the skeleton-based sweeping method[32,38] for mesh generation to represent the tree structures of neuron geometry. Given the geometry information of neurons, our method automatically reconstructs 2D network geometry with high accuracy and high order of continuity for IGA computation. Our automatic IGA optimization solver provides an efficient computation tool for studies of material transport regulation in complex neurite networks. The current 2D solver can be easily generalized to 3D and it is also extensible to solve other PDE-CO models of cellular processes in complex neurite network geometry.

Our study has its limitations, which we are addressing in the ongoing work. In the current model, we only consider the influence of traffic jams on the material concentration but neglect its effect on the deformation of neuron geometries. In addition, although IGA offers great advantages in accurately simulating material transport control in complex neuron geometries, the computational cost of simulating transport in large-scale neurite networks remains very expensive, which limits its biomedical application. To improve the computational efficiency of our model, we will adopt deep learning techniques to build fast and accurate surrogate models[50,73]. Moreover, designing comparable experiments is necessary to verify the active regulation mechanism in our model. For instance, the photoactivation technique[74,75] can be used to visualize the material transport process and extract the velocity or concentration distribution to compare with simulation results. The velocity data from experiments can also be used as the pre-defined velocities to provide a more realistic simulation. Regarding the control variables in our model, they may have relationship with the distribution of molecular motors that affect the transport dynamics on MTs, and thus needs further experimental data for verification. Despite these limitations, our simulation directly shows how the traffic jam is formed in neurons and how neurons could control material traffic to avoid traffic jams. The simulation results provide references to further answer the question of how neurons deliver the right material to the right destination in a balanced manner in their complex neurite networks and how the transport may be affected by disease conditions.

Supplementary Information.