Axonal transport during injury on a theoretical axon

1. Introduction

Traumatic brain injury (TBI) and the associated long-term disease progression known as chronic traumatic encephalopathy (CTE) continue to be the leading causes of mortality and morbidity worldwide, with implications for Alzheimer's disease and Parkinson's disease (Daneshvar et al., 2015). Trauma pathophysiology includes damage to neurons such as shearing, crush damage, and stretch injuries that result in axonal pathologies at the earliest stages of the disease. Distortion of axonal cytoskeleton structure and varicosities, which are axonal protrusions similar to beads on a string, may also be evident along axon length (McKee et al., 2013). Axonal transport is required for neuron homeostasis, and neural signaling and disruptions in the process are increasingly associated with a variety of neurological disorders (Sleigh et al., 2019). Shared neuropathological hallmarks such as the transactivation response (TAR) DNA-binding protein 43 (TDP-43) are seen not only in TBI but also in other neurodegenerative pathologies, such as Amyotrophic lateral sclerosis (ALS), for which axonal transport defects are common (Baldwin et al., 2016). TDP-43 is a multifunctional RNA-binding protein with roles in mRNA transport in axons and dendrites and contributes to localized translation in pre-synaptic neurons for rapid signaling and neuronal homeostasis (Ling, 2018; Nagano et al., 2020). Regulation of TDP-43 is linked to axonal survival and neurodegeneration in visual circuits (Shigeoka et al., 2016) as well as additional roles in the structure and function of dendritic spines in Drosophila and mice (Lu and Vogel, 2009; Fogarty et al., 2016; Handley et al., 2016). While disease hallmarks such as TDP-43 highlight the outcomes of disrupted axonal transport, theoretical models such as DyNAMO can now provide a detailed mechanistic understanding of the nano- to microscale processes that escape detection due to the limitations of live cell imaging at this scale. The DyNAMO theoretical axon is the first computational platform designed to enable studies of altered dynamics of multiple motor proteins in axonal injury scenarios, including visualization, classification, and quantification of transiting events and outputs to drive long-term advancements in treatments for TBI and other axonopathies.

In the three decades since the discovery of neuronal kinesin and fast axonal transport (Smith, 1980; Tsukita and Ishikawa, 1980; Vale et al., 1985), progress toward nanoscale visualization of dynamics along axons remains challenging (Surana et al., 2020). Biological studies of fast axonal transport in living cells, which initially relied on techniques such as video-enhanced contrast-differential interference contrast microscopy (Allen et al., 1982; Song et al., 2016), have evolved to apply correlative live cell and super-resolution microscopy (Bálint et al., 2013), as well as an advanced and growing toolkit of technical approaches (Surana et al., 2020), but still remain challenging for interpreting and evaluating interactions of multiple types of biomolecules at nanoscale resolution. In contrast, detailed imaging of fixed samples has led to a dramatic new realization of the complex and layered biological infrastructure of the axon that includes three interacting cytoskeletons, containing discontinuous, staggered, polar MT bundles (Yamada et al., 1971), actin waves (Ruthel and Banker, 1999), actin rings (Xu et al., 2013; Vassilopoulos et al., 2019), actin hotspots and trails (Ganguly et al., 2015), and neurofilaments (Trojanowski et al., 1986; Nixon and Shea, 1992; reviewed by Leterrier et al., 2017; Papandréou and Leterrier, 2018; Hahn et al., 2019; Leterrier, 2021), as well as numerous interacting motor and non-motor binding proteins (Hirokawa and Takemura, 2005; Conde and Cáceres, 2009; Hirokawa et al., 2009) essential to axonal functions. Advanced mathematical models are critical for bridging experimental gaps, such as nanoscale dynamics, while incorporating static information from high-resolution images on fixed samples or obtained by super-resolution methods on living cells into the theoretical heuristic. Few computational models exist to address axonopathies. A model allowing simulations of the viscoelasticity of the MT-associated protein tau, MAPT, was developed relevant to Alzheimer's disease and axonal stretch injuries (Ahmadzadeh et al., 2014, 2015; Barbier et al., 2019). Banerjee et al. (2018) modeled multiple neuron compartmentalized functions as a nanoscale end-to-end communication system, enabling experimental data from multiple biological studies to be evaluated to compare amyloid-beta neurotoxicity at membranes, receptors, and internal calcium stores relevant to Alzheimer's disease. The DyNAMO theoretical axon is a significant step in providing a detailed framework that significantly expands the ability to interrogate a diverse set of parameters relevant to axonal MT injury, infrastructure perturbations, and impacts on motor-based axonal transport relevant to neurological disease (Sleigh et al., 2019).

To model axonal transport, including non-equilibrium statistical mechanics and diffusion behaviors of multiple motors on axonal MT bundles, DyNAMO applies a totally asymmetric simple exclusion process, TASEP (Spohn, 1991; Derrida, 1998; Derrida and Evans, 2009), and Langmuir kinetics (LK). TASEP is a cross-disciplinary framework that restricts events to a unidirectional driven diffuse system relevant to diverse biological applications of intracellular molecular transport as well as the fields of chemistry and physics. LK algorithms, when integrated with TASEP, allow the capture of diffusive absorption/desorption dynamics, such as reattachment/detachment of motor proteins from MTs and inter-lane transition (Wang et al., 2007) as well as more complex scenarios in DyNAMO. Generalized transport studies of single or dual MTs have revealed the potential of combined TASEP-LK analysis on high- and low-density particle trafficking (Parmeggiani et al., 2003, 2004; Nishinari et al., 2005; Wang et al., 2007). DyNAMO significantly builds on these approaches to incorporate a rich set of parameters for dense capture of information, quantification, and co-evaluation relevant to axonal transport. This includes staggered discontinuous MT bundles, dynamic events of multiple motors with motor-specific motility parameters, and motor-type MT crowding, dissociation, and reassociation dynamics at site-specific locations along the theoretical axon, never co-evaluated. DyNAMO provides multiple types of quantifiable readings resulting from MT injury for the motor types navigating the theoretical axon, including input, stalling, progression patterns, and outputs. The roles and impact of varying kinesin types and their motility parameters in in vivo axon transport (Brunner et al., 2004; Verbrugge et al., 2009; Lo et al., 2011; Sun et al., 2011; Jenkins et al., 2012; Lessard et al., 2019), including axonal distribution (Hirokawa et al., 2009), remain unknown. Advanced theoretical axon models, such as DyNAMO, will benefit the understanding of axonal mechanisms and aid the analysis of brain injury pathologies and axonopathies that impact axonal MT structure and neuronal signaling. DyNAMO provides the ability to simulate MT length and multi-MT staggering injury changes in axonal MT cytoskeleton architecture along with impacts on the progression of multiple kinesins of differing motility and processivity. DyNAMO expands our understanding of axonal injury events and provides a means to spatio-temporally characterize and quantify associated transport impacts. This information and the continued advancement of neuronal theoretical models are anticipated to have numerous benefits, including mechanistic and biomarker discovery to distinguish minor vs. severe changes in axonal transport and other mechanistic considerations in designing future diagnostics and therapeutics for axonopathies.

2. Materials and methods

2.1. Development of a framework to capture complex axonal parameters

We developed a multiple MT lattice structure (MMLS) for axonal MT bundles that are designed on the basic framework of TASEP-LK, significantly modified to enable a nanoscale functional analysis of discontinuous staggered MTs with interacting dual kinesins and distinct kinesin-specific motility parameters in a theoretical axonal framework. Table 1 summarizes the experimental and defined parameters applied. Combined TASEP-LK has been applied by others to analyze unidirectional single motor flow in single MTs and motor particles for the co-existence of high- and low-density trafficking regions (Parmeggiani et al., 2003, 2004; Nishinari et al., 2005) or in two-MT models with motors as particles with symmetric or asymmetric distribution (Wang et al., 2007). In DyNAMO, the influx of the two kinesins into the MMLS at MT minus ends is regulatable and followed by stochastic progression of the kinesins along binding sites with lateral MT transitions permissible depending on MMLS scenarios. The only interactions considered between the two kinesins modeled in this study are those that result on a given single lane as hindered forward progression due to varied motility. To incorporate asymmetric coupling with MTs (Pronina and Kolomeisky, 2004, 2006; Jiang et al., 2007, 2008; Xiao et al., 2010), we used a three-MT MMLS and lateral movement parameters introduced for both kinesins to regulate crowding dynamics (Figure 1C). Two neuronal kinesins, Kin1 (S) and Kin3 (F) are described in the text, chosen for their significance and differences in motility and processivity (Verbrugge et al., 2009; Sun et al., 2011; Soppina et al., 2014; Lessard et al., 2019). LK detachment/reattachment dynamics (Jiang et al., 2007; Wang et al., 2007; Vuijk et al., 2015) are evaluated relative to kinesin-specific parameters and MMLS interactions. Detachment of kinesins into finite-length productive reservoirs at MT corresponding lattice sites in DyNAMO, instead of a singular bulk reservoir used in other studies, allows detailed retention of spatio-temporal dynamic information. The influx of kinesins into the MMLS is done via TASEP as bulk access to MTs and has measurable impacts on crowding and forward progression in coupled multilane interactions. Our DyNAMO model is the first theoretical axon model capable of capturing the nanoscale dynamics of multiple kinesins trafficking on a three-MT expandable TASEP-LK framework of discontinuous staggered MTs. DyNAMO expands the critical network of heuristic and mathematical analysis of compartmentalized sub-neuronal functions.

2.2. Numerical model of the axonal MT architecture

To represent the complex multiple MT and staggered architecture in the axon we extend a simple TASEP single MT framework, applying basic MT descriptive principles of a lattice with sites. According to Parmeggiani et al. (2004), a one-dimensional TASEP was used, in which the analytical solution applied mean-field approximation to the Heisenberg equation and compared outcomes with Monte Carlo simulation results. Similarly, we defined a finite one-dimensional lattice with sites labeled i = 1, 2, … N with lattice spacing l (8 nm). The stochastic equation for the bulk lattice sites (1 ≤ i ≤ N−1) with attachment (A) and detachment (D) dynamics of an MT is described by a master equation, where $\frac{d\langle n_l\rangle }{dt}$determines the rate of change in expected occupancy of a site and 〈 〉 characterizes the statistical average.

The subsequent work by Verma et al. (2015) and Verma and Gupta (2018) resulted in a heuristic designed for analysis of a single MT with two lanes. In DyNAMO the tracks represent MT protofilaments on MTs in fully asymmetric coupled TASEP-LK under open boundary conditions. The evolution of motor density in the bulk for two lanes was given by:

The positive and negative terms on the right-hand side denote gain and loss terms arising due to kinesin attachment, detachment, and forward hopping succeeded by lane-changing processes. At the boundaries, particle densities are given by:

2.3. Capturing individual motor variables in a combined system of multiple motors

To capture kinesin motility parameters along with any inter-dependencies of transport and processivity, we combined multiple parameters as step-dependent equations. The dependency on multiple motor parameters creates a challenge to incorporate all variables into solvable partial differential equations (PDEs). To be able to capture individual motor variables in a combined system of multiple motors, as well as any interdependencies, we combine multiple motility parameters of speed and processivity into the following step-dependent equations. We let $\langle n_t^i\rangle$and $\langle m_t^i\rangle$ denote the discrete binary states of occupancy of protofilament lattice sites by motor species “a” and “b” representing kinesins at site “i” in lane t, respectively. The forward hopping rate or step size of species “a” is normalized to be a unit, and the forward hopping rate or step size of “b” is assumed v ≤ 1. The reattachment (R) rate of species “a” and “b” to any unoccupied lattice site i is denoted by ω_{a, R} and ω_{b, R}, and the detachment (D) rate from the lattice site is denoted by ω_{a, D} and ω_{b, D}. For sites 1 ≤ i ≤ N−1, the time evolution of $n_t^{i}and{m}_t^i$are governed by the following PDEs to incorporate asymmetric coupling of motors in the lanes with TASEP-LK and limited processive steps on lane t:

The complex correlation between multiple motors in the MMLS model made it impossible to find an analytical or numerical solution using partial differential equations (PDEs) without removing these interdependencies. Furthermore, external parametric conditions such as the processivity and lifetime of motors added more complexity to the equation. Thus, a heuristic approach was used to solve MMLS, which involved defining specific testable algorithm conditions to examine these complex interactions. The forward hopping of motors and asymmetric coupling with adjacent lanes is governed by the TASEP principle but our model allowed for motors to make multiple hops with their respective speeds within one-time frame. The natural exit and entry of motors from/into the productive reservoir into/from the MT lanes are governed by the principles of LK. The entire dynamics of the motor are subdivided into four major action phases—loading, forward progression, detachment–reattachment, and output/delivery governed by the TASEP-LK principles with asymmetric coupling.

2.4. The influx of kinesins onto the minus ends of MT MMLS lanes

In DyNAMO, the finite junction reservoir loads Kin3 (F) and Kin1 (S) onto MT binding sites following a normalized Gaussian distribution. The inter-arrival time or the time between successive loading is distributed. The expected run time of the simulation is proportional to motor generation rate α_a, α_b (“influx rate” in motors/s). In Scenarios 1–8, both Kin3 (F) and Kin1 (S) are loaded onto MTs one kinesin at a time from the finite productivity reservoir present at the minus end of the MT lanes. As described in Algorithm 1, kinesin loading is restricted to the binding site i = 1 (referred to as the loading site) at the minus end of each MT. The kinesins have access to all parallel lanes with equal probability in all scenarios, except when the flow is forcefully channelized to a single lane due to proximal staggering in Scenarios 3–5.

2.5. Unidirectional traversal of kinesins along the MT lane (TASEP)

We confined our MMLS to a two-dimensional plane, which is necessary to simplify the current simulation to implement greater overall complexity with additional parameters and realize that this may create some bias (Whitesides, 1970). In our two-kinesin DyNAMO trafficking model, we have defined multiple one-dimensional LS of T lanes (T ≥ 2) of length L consisting of N lattice sites (N ≥ 3). The multiple one-dimensional binding sites are labeled i = 1, 2, 3 ... N, and lattice spacing is l=L/N. The binding sites i = 1 and i = N define the left and right boundary, while i =2, 3 … N-1 is referred to as the bulk. The binding site is characterized by $\left\{n_{t=1,2,3\dots t}^{i=1,2,3\dots N}\right\}$, when the occupied state by motor “a” (denoted by $n_t^i$) or motor “b” (denoted by $m_t^i$) is equal to zero (vacant) or one (occupied). Every site can be occupied by either kinesin or not occupied by either kinesin type.

To model the transport of two different types of motors, the faster Kin3 (F) can traverse to the next open site at a rate of v_a steps, while the slower Kin1 (S) can traverse at a rate of v_b steps (assuming: v_a ≥ v_b). No kinesin can hop over another to occupy a forward available site during traversing along the lane (protofilament). Kin3 (F) upon encountering a barrier (on the occupancy of the next site by Kin1 (S) or blockage), either laterally moves into an available lateral MT lane or gets detached into the productive reservoir. Kin1 (S) is a limited processivity motor. To characterize the processivity component of motility, we implement Kin1 (S) with a processivity length (the maximum distance the motor traverses on an MT before detachment) of 1,070 nm (i.e., equivalent to ~130 lattice sites). The progressive system dynamics for different MT transitions based on the modified TASEP are governed by the following methodologies:

1. Loading: At the site (i = 1), a motor can enter the MMLS with rate α if unoccupied.

2. Forward progression: If the sites (i+1, i+2 … i+v) are unoccupied, a motor with velocity v can jump from site (i) to site (i+v) that is equivalent to an 8 nm dimeric kinesin step in DyNAMO.

3. Lateral movement: If any of the sites (i+1, i+2 … i+v) are occupied, a motor can jump from site (i) to either site (i-1) or (i+1) of an adjacent lane if unoccupied.

4. Delivery/output: At the site (i = N-v+1, N-v+2 … N), if a motor with velocity v is present it can leave the MMLS from the plus end of the MT lane.

Limited processivity requires the movement of kinesins to detach from MTs into an inter-MT reservoir from which they can either reattach onto the same MT or an adjacent MT in the MMLS. When congestion prevents reattachment to the same MT, reattachment must occur to an adjacent MT or the kinesin is retained in the productive reservoir until sites are available. A general Algorithm 2 governing the TASEP model is formulated on the basis of the above methodologies. Different simulation paradigms are considered for different scenarios. For example, in Scenario 2, we limit the provisions for lateral movement; whereas in all other scenarios, both Kin3 (F) and Kin1 (S) have provisions to move to adjacent MTs upon encountering congestion. In Scenarios 3–8, a portion of MMLS is staggered to channel kinesin flow into a singular MT lane either in proximal or distal ends to evaluate blocked access by a variety of mechanisms (damage, tubulin modifications, MAPs, and so forth) as described in the text. The kinesins can overcome crowding by detaching and reattaching to the localized productive reservoir and MT availability for movement to an adjacent lateral MT. Previous studies have included similar concepts of MT lanes and lattice sites (Dixit et al., 2008; Che et al., 2016; Liang et al., 2016), productive reservoirs (Leduc et al., 2012; Ciandrini et al., 2014; Feng et al., 2018), or side-stepping lateral movement to an adjacent MT lane (Wang et al., 2007; Jiang et al., 2008; Hoeprich et al., 2014), which are also conceptually included in DyNAMO.

2.6. Reattachment and detachment dynamics into productive reservoirs (LK)

Productive reservoirs are defined by LK from the lattice sites (1 ≤ i ≤ N-1) of the MMLS. In previous studies, TASEP-LK (Wang et al., 2007; Vuijk et al., 2015; Dhiman and Gupta, 2016) was used to define the reattachment and detachment of kinesins to and from a bulk productive reservoir. In DyNAMO, we define productive reservoirs instead as a queue of a fixed size corresponding to each lattice site (motor binding site) on MTs of the MMLS. The LK dynamics of reattachment and detachment are defined accordingly:

1. Detachment: During the kinesin motility along one-dimensional lattice sites, detachment occurs when one motor cannot take the assigned hops (stepping) due to crowding. That is, if the site (i+1, i+2 … i+v) of the MT lanes, such as Kin3 (F) trailing Kin1 (S), will result in Kin3 (F) being detached into a productive reservoir at lattice site i.

2. Reattachment: In the productive reservoir at lattice site i, within a given lifetime of the motor, if a site (i-1, i, i+1) is unoccupied, then the motor reattaches to an open site on the MT.

3. Leakage: If a motor does not reattach within its lifetime, then we describe the motor as leaked out of the system into a non-productive reservoir.

Algorithm 3 describes the process of reattachment/detachment of motors (also described as absorption and desorption). We apply LK to describe the off–on dynamics of motors along the MTs because of crowding, which also creates opportunities for motors to reattach to lateral MTs to avoid local crowding along MTs. All our scenarios considered for simulation have provision for kinesin reattachment back to MTs on the availability of vacant sites corresponding to the point of detachment. Scenarios 1 and 2 featured parallel MT lanes, while Scenarios 3–8 included staggering to limit access to adjacent MTs, either proximally at the minus end or distally at the plus end of MTs. The LK model enabled the interpretation of the on–off dynamics of kinesins under different scenarios and flow rates.

2.7. Delivery of kinesins from the MT plus end

The transit of Kin3 (F) and Kin1 (S) is to the distal plus end (i = N) of an MMLS by motor service rate β_a, β_b given as an “Outflow rate” in motors/s). At the end of an MT lane, the kinesins are transferred to a junction reservoir (described in Algorithm 4). These kinesins at the junction reservoir form a limited pool of available kinesins that can be loaded onto the next MMLS segment. This limited reservoir of kinesins is analyzed in the evaluations of mixed organization of MMLS scenarios in an extended axonal framework.

2.8. Solution scheme paradigms

The simulation environment uses Python 3.8 for testing the algorithms in DyNAMO. The code and documentation developed for DyNAMO are available in the open platform GitHub. The parameters used are in accordance with biologically observed data obtained from numerous different sources and summarized in Table 1. In summary, the model starts by setting up the MMLS system of T number of parallel lanes and N binding sites. The dual-motor dynamicity is parameterized by their individual speed/velocity (v), processive steps (p), and influx rate (α). Based on the chosen scenario, the MMLS is set up with different boundary conditions and restrictions within the system. The kinesins from the finite reservoir in the minus end are loaded onto MT initial sites at each iterative time stamp. The motors based on the different decisions governed by TASEP and LK, gradually traverse along MTs and exit the system from MT plus ends. The evaluated metrics (such as outflow rate, dissociation, and dynamicity) and critical condition cases are considered for depicting different kinesins of distinct parameters trafficking across various scenarios. To extensively study the effect of these parameters on multiple linked segments of MMLS scenarios, we considered a linked architecture with varied linked MMLS having different scenario restrictions. These individual MMLS segments evaluated in a multiple MMLS scenario retain their outcomes as discrete segments but when linked demonstrate the broader impact of distributed or concentrated MT injuries. The simulations were run on standard desktop processors (Intel i7-9th generation with RTX 2080 Ti) and took ~3–4 h for each scenario. The subcases of each scenario were also evaluated on GPU clusters at Purdue University. In future, the DyNAMO model can be scaled up to analyze longer MT tracks with more protofilaments per MT, as well as incorporate additional kinesins. However, such an increase in variables will significantly increase the computational complexity of the model and require access to more powerful computing resources such as larger clusters or supercomputers.

3. Results

3.1. DyNAMO, a parameter-dense theoretical axon framework

To evaluate combined contributions of multiple kinesins on a disease- or injury- perturbed axonal MT infrastructure requires the ability to capture detail-rich dynamic information. In DyNAMO, we apply a TASEP-LK paradigm approach for non-equilibrium and equilibrium statistical mechanics, summarized here and detailed in Methods. The simulations are analyzed in the environment of Python 3.8, with parametric values as defined from experimental measurements or by assumptions based on published experimental data (Table 1). Two unidirectionally trafficking kinesins with different motility and processivity rates are mathematically described to enable considerations of crowding, including influx rates, and the influence of motor type in damage scenarios. In this study, we focus on human conventional Kinesin-1 KHC and Kinesin-3 Kif1A neuronal kinesins referred to herein as Kin1 (S) and Kin3 (F) slower and faster motors. Electron microscopy studies of longitudinal or transverse axon sections were used to provide guidance on parameters for MT discontinuous organization. In human hippocampal developing neurons, shorter MT lengths are reported (4.02 +/– 5.28 um; Yu and Baas, 1994). Detailed MT cytoskeleton maps along axon length in Caenorhabditis elegans neurons are generally slightly longer 1.2–10.7 um in length (Chalfie and Thomson, 1979). We apply a three-MT framework MMLS during simulations with the ability of kinesins to move between all MTs and capture features of normal or increased staggering as a simulation of MT injury (Matamoros and Baas, 2016; Figure 1). The MT surface is modeled as a single lattice structure, single protofilament-like, with N kinesin-binding sites along a defined but adjustable MT length. Axonal MTs are non-contiguous staggered short segments (Figure 1A). Maximum MT lengths for this study are set at 4,000 nm, reflective of average MT segment lengths observed experimentally (Yu and Baas, 1994), and containing up to 500 binding sites of stochastically Gaussian-envelope-like kinesin steps of 8 nm, the step size of dimeric kinesin moving along a protofilament (Svoboda et al., 1993; Kojima et al., 1997; Coy et al., 1999; Okada et al., 2003).

In TASEP, motor interactions with MTs are described by hard-core repulsion, such that each binding site is accessible to only one of any kinesin at any moment. TASEP describes bulk kinesins entering or leaving the system at system boundaries (MT ends), such as from the AIS at the neuron soma (Petersen et al., 2014; Kapitein and Hoogenraad, 2015; Leterrier, 2016) or in interconnected multiple MMLS along the axon length. LK allows specified actions on transiting kinesins to be mathematically described, such as kinesins undergoing MT-to-MT jumps or dissociating into and reassociation from position-correlated inter-MT reservoirs. Detachments can occur via internal crowding, processivity limitations, or termination of MT segments, with MT reassociation when there is open access. The use of positionally distinct inter-MT reservoirs in DyNAMO rather than a long bulk reservoir (Leduc et al., 2012; Feng et al., 2018) allows us to capture stationary positional with on/off dynamicity. High dynamicity reflects crowding and highlights delayed forward progression points along the theoretical axon. The positional reservoirs have an adjustable storage capacity that can be restricted in capacity and time. Kin1 (S) has been observed to have a long but definitive lifespan of ~37 h in the vicinity of axonal MTs (Brunner et al., 2004; Verbrugge et al., 2009), and in our shorter-term scenarios, no motors are lost from inter-MT reservoirs, and all eventually return to the system. The full simulation environment is detailed in Methods. The corresponding kinetic rates obey a detailed balance such that the system evolves into an equilibrium steady state.

3.2. Initial patterning and progression of kinesins on the theoretical axon

In previous studies applying TASEP to two-lane systems, single particle dynamics vary due to exchange between lanes at asymmetric rates (Pronina and Kolomeisky, 2006; Jiang et al., 2007, 2008; Cook et al., 2009; Vuijk et al., 2015; Dhiman and Gupta, 2016; Gupta, 2016). Using DyNAMO, we evaluated the impact of “influx rate” in our multi-kinesin MMLS to study how kinesin-specific motility parameters such as speed and processivity, as well as MT-to-MT lateral access or restriction, impact progression. Two MMLS scenarios and outcomes are described in Figure 2. In MMLS Scenario 1 (Figure 2A), all MTs are accessible with lateral movement at any position. This allowance accommodates the presence of Kin1 (S) and Kin3 (F) of variable speed by allowing Kin3 (F) to jump laterally to adjacent MTs to avoid crowding and support forward progression. In Scenario 2 (Figure 2B), no lateral access is allowed. When the Influx rate is varied incrementally from 10 to 15 motors/s, Kin1 (S) is minimally affected, whereas the impact on the detachment of Kin3 (F) is significant with 400–700 motors detached per second. Comparative histograms (Figure 2C vs. Figure 2D) and quantitative analysis (Figure 2E vs. Figure 2F) are shown. The primary contributing factor to the detachment of Kin3 (F) is crowding, whereas, for Kin1 (S), it is limited processivity, which is consistent with previous theoretical models (Leduc et al., 2012; Ciandrini et al., 2014). In our analysis, early lateral access to MTs in axons has the primary impact on defining the initial kinesin traffic patterning in a multiple motor system, which is consistent with experimental observations in vivo (Rank and Frey, 2018).

3.3. Simulated axonal injury differently impacts kinesin types in dissociation or progression

To simulate axonal injury such as stretch-induced gaps or destabilization of MTs affecting length due to partial depolymerization, we modeled the impact of altered MT staggered lengths in our MMLS on kinesin progression vs. on/off dynamicity. Proximal or distal staggering was evaluated in multiple MMLS scenarios with 12.5% staggering considered within normal variation (Scenario 3) and 25–50% staggering reflecting increasing degrees of injury-induced perturbation (Scenarios 4 and 5, respectively). Staggering creates channelization of kinesins (Figures 3A–C). When staggering is proximal and the influx rate is low (< 13 motors/s), Kin3 (F) is impacted, whereas Kin1 (S) is not impacted (Figures 3D, E). However, at higher influx rates (>13 motors/s), bulk crowding occurs due to the channelization of multiple kinesin types with differing motility parameters onto a single MT lane. This results in a high rate of detachment of both Kin1 (S) and Kin3 (F) into the inter-MT reservoir. The detachment rate for staggered injured MT tracks (such as Scenarios 4 and 5) increases by 3–4 times compared to a normal MT Scenario 3 for faster Kin3 (F), while it results in an almost 10-fold increase for the slower Kin1(S). In addition, channelization and influx rates combine to limit the availability of free MT sites for kinesins to move along a singular MT lane. Thus, a huge dynamic on-off activity is observed with an increased length of channelization. DyNAMO reveals that during proximal channelization, and dependent on influx rate, the slowest moving kinesin dictates the progression of all the kinesins (see Figures 3F, G) until additional parallel MT lanes become available for MT-to-MT jumps. The outcomes support a mechanistic need to limit staggering within the AIS or apply gating controls at the AIS that have been proposed biologically (Petersen et al., 2014; Kapitein and Hoogenraad, 2015; Leterrier, 2016).

Increased staggering length at the distal end of the MMLS (Figures 4A–C) was also modeled at 12.5% (Scenario 6), 25% (Scenario 7), and 50% (Scenario 8). In distal staggering, the initial uniform loading and distribution of the motors along the parallel MT lanes abruptly transitions to channeling. Kin1 (S) and Kin3 (F) motors become channelized onto a dedicated single MT lane. At a low influx rate, Kin3 (F) is primarily affected, whereas a higher density of motors results in bulk crowding that impacts both Kin1 (S) and Kin3 (F) forcing increased detachment. Histogram plots of staggered vs. open MT regions (Figures 4D, E) show the effect of distal channelization on kinesins with different motility rates. In summary, an increased initial influx rate of motors from 10 to 15 motors/s results in denser packing of motors on the MMLS that increases dynamicity abruptly at the junctional change to staggering. The channelization length (12.5, 25, or 50% blockage of lateral binding sites) impacts the progression rate, defined by the slowest motor, and further impacts the degree of dynamicity along the staggered column (Figures 4F, G). The rate of on–off activity increases proportionately with increases in the length of staggered columns. The DyNAMO model reveals the importance of mapping the natural and injury distribution of staggered axonal MTs along the axon length in future. Furthermore, DyNAMO reveals that proximal staggering of axonal MTs at the AIS or along the axonal MT length can generate a gatekeeper effect dominated by slower kinesins, whereas distal staggering has a clearing effect that can remove the bulk of kinesins such that their reattachment will immediately influence the downstream distribution of kinesins on the theoretical axon and therefore signaling output.

3.4. Kinesin dynamicity is elevated in axonal injury and impacts concentrations of axon-distributed motors

We refer to the detachment and reassociation of kinesin as dynamicity and compare different scenarios of MT lattice staggering in the MMLS and influx rates on the dynamicity of kinesin types (Figure 5). The DyNAMO allows kinesin reattachment from reservoirs to MTs at plus/minus one step around that detachment site. Four MMLS scenarios with influx rates of 10 or 15 motors/s were compared (Figures 5A–D), namely Scenario 1 (maximum MT lengths and MT-to-MT lateral access), Scenario 2 (maximum MT lengths with no lateral movement permitted), Scenario 4 (25% proximal staggering), and Scenario 7 (25% distal staggering). For each scenario and influx rate, four dynamicity maps are shown that represent data for multiple conditions that are Kin3 (F) and Kin1 (S) separately, both motors combined, and Kin1 (S) alone on an expanded x-axis scale. In Scenario 1, which is least restrictive for MT access, the rate of detachment for both Kin1 (S) and Kin3 (F) is nominal at 10 motors/s (Figure 5A, central column, note scales, and dispersed peaks). The Kin3 (F) processivity is equivalent to the MT maximum length in our defined parameters, and the allowance of MT-to-MT transitions results in minimal dissociations affected by Kin1 (S) in Scenario 1 at 10 motors/s. At an increased 15 motor/s, the dynamicity of both Kin1 (S) and Kin3 (F) increases due to a higher crowding of motors in between MT lanes (Figure 5A, right, note scales, and compacting peaks). In Scenario 2, when no lateral MT-to-MT access is permitted (Figure 5B), the flow rate and composition of bound motors on MTs are dictated by slower Kin1 (S). This is due to the forced dissociation of Kin3 (F) when it trails and encounters Kin1 (S) (Figure 5B central column, note scales, and elevated compacted dynamicity plot). The elevated blue peaks are in accordance with the high dynamic on-off characteristic of Kin3 (F) along the bulk. In scenarios with proximal (Scenario 4, Figure 5C) or distal staggering (Scenario 7, Figure 5D) and a slower flow rate of 10 motors/s, the Kin1 (S) motor is nominally affected, whereas Kin3 (F) due to hindrance along the singular MT lane by the slower Kin3 (S) has increased dissociation to the productive reservoir and dynamic on–off flow in channeled regions. At 15 motors/s, both kinesins are limited in progression, with increased dissociation and limited ability to reassociate due to crowded conditions, effectively gating the downstream flow of motors and the type of motors progressing (Figure 5C, central and right columns, note scales). In situations involving distal staggering and high influx rates of 15 motors/s (Figure 5D, right, note scale), we observed bulk detachment of both motors, resulting in a much crowded and channelized flow with limited progressive movement due to limited access to rejoin MTs. This is reflected in a flatter dynamic plot, indicating a narrow band of flow. The distribution and flow dynamics of the kinesins along an MMLS and productive reservoir in two staggered scenarios (Scenarios 4 and 7) are shown in Supplementary Videos 1, 2. Overall, our findings suggest that under high flow rates of 15 motors/s, the channelization of motors at the distal end of the MTs affects the dissociation of both faster and slower motors. The results of our DyNAMO simulations highlight that the alignment of axonal MT organization with different kinesin parameters, such as speed and processivity, plays an important role in the spatial flow of kinesin types and concentrations along the axonal MT network.

3.5. DyNAMO analysis of multiple injury scenarios along the axon length and signal outputs

To gain a better understanding of how motors flow in bulk along the axon in different scenarios and segments, we developed a multi-MMLS heuristic consisting of three linked MMLS scenarios and a junctional reservoir. We computed the outflow rate of each motor from each individual lane/MT of the MMLS to assess the impact of different scenarios on the output flow of Kin3 (F) and Kin1 (S). We show two methods for representing motor delivery at different rates: “DETAILED” flow, which shows motor delivery at every computation timestep (40 ms; Figures 6A, B, upper rows), and “BINNED” data, which is a collective output over a timeframe of 1 s, equivalent to 25 iterations (Figures 6A, B, bottom rows). To illustrate the temporal output dynamics, we also generated simulation videos (Supplementary Videos 3, 4) for specific MMLS Scenarios 4 and 7. The temporal output of the two motor types was tracked through each of the single scenario MT sections of the multi-MMLS scenarios from left to right. Each dot along the time (x) axis represents the cumulative number of motors output from each MT section per second (values indicated on the y-axis; Figures 6A, B, right). In the multi-MMLS Scenario 1-4P-7D (P-25% proximal staggering and D-25% distal staggering), the binned plot in Section 1 reveals that motors are equally output from all three lanes due to the provision of lateral movement on adjacent lanes that avoids crowding. As motors enter and exit Section 2, proximal staggering exists, and the slower motor Kin1 (S) dominates the singular channeled lane (Lane 2), while faster motors tend to move to adjacent lanes (Lanes 1 and 3). The output of Section 2 reflects the pattern established by the proximal staggering, as motors retain their favorable positions. Entering and exiting Section 3, since all the outflows are channelized at the distal end to a single lane (Lane 2), the output reflects a restored mixed motor composition. When a different MMLS Scenario 1-7D-4P is evaluated in which distal and proximal staggering are adjacent, the final output from each MT is reflective of the previous injury event. In both evaluated multi-MMLS scenarios, we observe that the temporal outflow rate of kinesins is ~4–5 motors/s for all sections. To further investigate the relationship of kinesin processivity and MT length with the temporal outflow, we simulated similar multi-MMLS scenarios with a longer MT length of 8 μm (Supplementary Figure 1). We observe that for longer MT lengths within the multi-MMLS, the temporal outflow rate is comparatively lower (~2–3 motors/s), an indication that MT displacements of kinesins, on/off dynamicity, and loss of kinesins from the productive reservoir are greater.

In Figure 6, we evaluated the motor distribution output for different multi-MMLS scenarios under a continuous motor flow rate of 10 motors/s. We next evaluated the impact of a pulsed input of motors in several multi-MMLS scenarios at 10/15 motors/s. We either generate a “fresh” block and release it to free MTs (Figure 7B) or we assume MTs already contain residual motors “sprinkled” due to “previous” unblocked trafficking in that region (Figure 7C). The modeled distributed P-D multi-MMLS scenarios are 1-3P-6D (12.5%), 1-4P-7D (25%) (Figure 7A), and 1-5P-8D (50%), and D-P multi-MMLS scenarios that create a larger injury site are 1-6D-3P (12.5% staggering), 1-7D-4P (25% staggering), and 1-8D-6P (50% staggering). All injury scenarios (25 or 50% staggering) are compared to the 12.5% staggering considered within the normal range (Figures 7B, C, gray rows). We applied a pulsed input for 15 s which involved 300–450 motors for both Kin1 (S) and Kin3 (F) at association rates of 10 and 15 motors per second, respectively. We recorded the time it took for all the motors to exit Section 3 (pulsed time output) in each multi-MMLS scenario, and the results are listed in the tables (Figures 7B, C). The green-shaded ranges are near-normal (gray shading). When intermittent successive regions of staggering are present, such as in multi-MMLS P-D scenarios, kinesins take longer to exit (delayed progression) than when a larger single continuous staggered region (our multi-MMLS D-P scenarios) is present. Thus, DyNAMO predicts that smaller intermittent damage is more disruptive in injury scenarios. The higher flow rate in injury scenarios coupled with pulsed input in a sprinkled multi-MMLS, generates the greatest delays in output for both motors. Kin3 (F) is affected by “full” block and release for 15 motors/s, while Kin3 (F) and Kin1 (S) are also affected by slower inflow in MMLS for a higher percentage of staggering. All these data are consistent with counteracting factors of motor speed, processivity, MT organization, and crowding.

4. Discussion

Mathematical models of neuronal axons complement and expand biological studies by providing a missing correlation between nanoscale mechanisms and spatio-temporal dynamic events relevant to neuronal function and signaling. The nanoscale dynamics of motors along the axon's dense polar (Burton and Paige, 1981) naturally staggered and discontinuous MT networks (Baas et al., 2016; Prokop, 2020) in living biological systems are still largely undescribed experimentally. This limits the full understanding of the normal and perturbed dynamic mechanisms of multiple interacting biomolecules in living systems. Moreover, surprisingly few studies exist that have detailed regional maps of axonal MT organization (Chalfie and Thomson, 1979; Yu and Baas, 1994; Reis et al., 2012), and virtually no study exists that map axonal MT organization and staggering distribution along the entire length. Due to limited nanoscale experimental correlations that reflect dynamics, axonopathies and axonal transport deficits in neurological disorders remain largely mysterious, hindering our ability to unravel neuropathology mechanisms (Millecamps and Julien, 2013). The DyNAMO theoretical model of axonal transport was developed to reveal the complex interrelationship between normal and injury-perturbed MT axonal architecture and the impacts on interactions of multiple motor protein types traversing that infrastructure. DyNAMO permits multi-parameter co-evaluations relevant to axonal injury (Figure 7D), including multiple kinesin types with distinct speed and processivity motility parameters, the ability to influence influx rate, and to track motor movements as MT-to-MT, detachment/reassociation dynamicity, stationary states, or forward progression, and in the context of axonal MT injuries reflected in MT length (or blocked access) and MT staggering that are regional or distributed lengthwise. No other model currently exists with this complexity to simulate axonal transport.

DyNAMO provides modeling of axonal injury and is expandable to consider broader axonal mechanisms and applications. In the context of TBI and CTE, DyNAMO is expected to help identify, classify, and model disease hallmarks and biomarkers relevant to injury severity by kinesins and their cargos of vesicles, organelles, RNA, proteins, and neurotransmitter signaling molecules, For example, recent studies indicate that TDP-43, an indicator of severe damage to axons, is transported by KHC in neurons, Kin1 (S) modeled here, with the overexpression or inhibition of KHC leading to an increase or decrease, respectively, in the transport of TDP-43 (Chu et al., 2019). DyNAMO analysis is consistent with this finding and reveals an impressive tolerance of slower Kin1 (S) kinesin to minor alterations in MT infrastructure, with effects only seen with more severe damage scenarios, supporting TDP-43 as a severe injury hallmark. In contrast, we observe immediate effects with minor damage on faster kinesins such as Kin3 (F) and their cargos. In relevance to signaling output, DyNAMO reveals how axonal injuries can delay signaling and generate altered signaling patterns, composition, and strengths in relevance to kinesin types. DyNAMO allows modeling multiple kinesins in relation to perturbations to axon MT infrastructure in this study, relevant to injury states (Matamoros and Baas, 2016; Muñoz-Lasso et al., 2020) but is generalizable to address a multitude of future questions on additional mechanisms that impact kinesin access to MTs in the context of MT-associated proteins (MAPs), such as CAMSAP (Jiang et al., 2014; Yau et al., 2014) and MAPT (Barbier et al., 2019), or post-translational tubulin modifications (Janke and Chloë Bulinski, 2012) postulated to represent a code for binding interactions, as well as MT length regulating agents (Conze et al., 2022). We also demonstrate DyNAMO's flexibility to evaluate additional parameters, such as reduced local levels of ATP available to kinesins in injury zones, a consequence of disrupted mitochondrial transport in neurodegeneration (Sleigh et al., 2019; Supplementary Figure 2). By subclassifying our position-specific productive reservoir, we regionally evaluated ATP concentration differences relevant to injury and non-injury sites. In the case of low ATP concentration, we defined a longer time constraint on kinesins vs. no time restriction on normal ATP levels in non-injury regions. The simulation output reveals that when the availability of ATP is reduced, there is delayed output flow of both motor types, with a higher impact on Kin1 (S) low processivity kinesin, and with output further influenced by the damage scenarios. In DyNAMO, any impact that removes a motor from the MT to the productive reservoir (i.e., crowding or limit of processivity or lack of ATP) offers an opportunity to control motor stationary position vs. forward motion, or in future studies, a switch between the anterograde and retrograde movement of cargo. Although in this study we focus on kinesin anterograde transport and MT axonal injury scenarios, intracellular transport of cargo is often bidirectional. This is due in part to retaining a mixed steady-state population of opposite polarity unidirectional kinesin and dynein motors (Encalada et al., 2011). Simulations will help distinguish multiple mechanisms that may include the ratio of motor types in a tug of war (Gross, 2004; Welte, 2004) or regulation of motor activity, including autoinhibition (Verhey and Hammond, 2009; Akhmanova and Hammer, 2010) or KLC interactions between the intact Kinesin-1C/KLC1 complex and dynein, which were shown to be needed for proper dynein retrograde activity (Martin et al., 1999; Ally et al., 2009; Encalada et al., 2011). Since a switch to the retrograde motion of motor-coated vesicles may impact several local levels of kinesins, consideration of the influx of the new kinesins will also need to be evaluated. In our studies with DyNAMO, we note that the rate of influx of motors has interestingly different impacts on kinesin types. We also anticipate the use of DyNAMO in the analysis of additional neuronal kinesins (Hirokawa and Takemura, 2005; Hirokawa et al., 2009) or evaluating the impact of multiple kinesins on processivity and obstacle avoidance (Ferro et al., 2019). DyNAMO reveals how axonal MT organization and multiple kinesins interplay to generate localized axonal crowding and patterning. Actin rings have also been proposed to play a role in localized crowding along the axon length (Sood et al., 2018). To what degree these two mechanisms, axonal MT staggering and actin rings, may co-contribute to the axonal flow of information is not yet known. In summary, DyNAMO is a powerful expandable platform for rigorous investigation of axonal transport mechanisms that reflect biological mechanisms and disease.

Finally, we anticipate that DyNAMO and other sub-neuronal models will also benefit the composition of larger neuronal communication models. We previously modeled the neuron as a nanoscale communication system using experimental data from multiple biological studies to understand and compare amyloid-beta toxicity and altered calcium signaling in Alzheimer's disease (Banerjee et al., 2018) in end-to-end events in intracellular signaling. It is anticipated that neuronal compartment models such as DyNAMO for axonopathies will advance biomedical research, analogous to models of synapses to evaluate synaptopathies for autism spectrum disorders (Chatterjee et al., 2021) and the spatial dendritic context of synaptic signaling (Larkum and Nevian, 2008; Tønnesen and Nägerl, 2016) in cognition and psychiatry. New initiatives toward mapping the spatial three-dimensional cellular context of rodents and eventually the human brain also make it imperative to continue to build and strengthen models of neuronal function (McMenamin et al., 2003; Watanabe et al., 2004; Pa?ca, 2018). Fully developed software and hardware computational models that interlink synaptic function, dendritic mechanisms, somatic correlations, and axonal signaling are almost certain to be as complex as larger topological models of brain function (Hennig, 2013; Chen and Sneyd, 2015). DyNAMO represents an innovative significant step in deciphering, visualizing, and quantifying axonal MT cytoskeleton end-to-end transport mechanisms in order to keep pace with the experimental realization of nanoscale structures and processes in biological systems. DyNAMO use will benefit studies relevant to a range of axonopathy diseases and mechanisms, including relevance to brain trauma or cognitive disorders such as Alzheimer's disease. The greatest challenge in modeling biological aspects of axonal transport with DyNAMO is expected to primarily be in computing power; currently, simulations are run on a laptop or small GPU cluster, whereas more complex parameters and interrelationships are expected to expand this need, particularly for full three-dimensional complex interrelationships. Where experimental detail is lacking, assumptions must be generated that may have compounding downstream effects. Because our focus was on the interrelationship of MT infrastructure and kinesin types, layers of kinesin regulation were not yet modeled in DyNAMO, such as multimers of kinesins that can affect processivity and bypassing obstacles, as well as dimerization that affects processivity. DyNAMO establishes a robust platform that can be built on by the user to make parameter and simulation decisions based on mechanistic questions.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material. The DyNAMO algorithms can be found in online repositories at this link: https://github.com/soumyadeepchandra/DyNAMO.

Author contributions

JP and AM conceived the project's experimental design. SC was primarily responsible for software algorithms and scenario data analysis. SC, ZO, and RC composed figures. JP, AM, and SC wrote most of the manuscript with contributions from all authors. All authors provided input to the mathematical framework, figures, and tables. All authors contributed to the article and approved the submitted version.

Funding

The Empire State Development's Division of Science, Technology, and Innovation (NYSTAR) Division funds 15 Centers for Advanced Technology (CATs) to encourage greater collaboration between private industry and universities in the development and application of new technologies. The CAT program, created in 1983, facilitates a continuing program of basic and applied research, development and technology transfer in multiple technological areas, in collaboration with and through the support of private industry. CATs play a critical role in spurring technology-based applied research and economic development in the state; promoting national and international research collaboration and innovation; and leveraging New York's research expertise and funding with investments from the federal government, foundations, businesses, venture capital firms and other entities. The work was supported by Center for Advanced Technology in Nanotechnology and Nanoelectronics funds to JP (CATN2 1186199-1-92476) as well as multi-university resources in a collaborative effort.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncel.2023.1215945/full#supplementary-material

References