Physical Biology of Axonal Damage

Excessive physical impacts to the head have direct implications on the structural integrity at the axonal level. Increasing evidence suggests that tau, an intrinsically disordered protein that stabilizes axonal microtubules, plays a critical role in the physical biology of axonal injury. However, the precise mechanisms of axonal damage remain incompletely understood. Here we propose a biophysical model of the axon to correlate the dynamic behavior of individual tau proteins under external physical forces to the evolution of axonal damage. To propagate damage across the scales, we adopt a consistent three-step strategy: First, we characterize the axonal response to external stretches and stretch rates for varying tau crosslink bond strengths using a discrete axonal damage model. Then, for each combination of stretch rates and bond strengths, we average the axonal force-stretch response of n = 10 discrete simulations, from which we derive and calibrate a homogenized constitutive model. Finally, we embed this homogenized model into a continuum axonal damage model of [1-d]-type in which d is a scalar damage parameter that is driven by the axonal stretch and stretch rate. We demonstrate that axonal damage emerges naturally from the interplay of physical forces and biological crosslinking. Our study reveals an emergent feature of the crosslink dynamics: With increasing loading rate, the axonal failure stretch increases, but axonal damage evolves earlier in time. For a wide range of physical stretch rates, from 0.1 to 10 /s, and biological bond strengths, from 1 to 100 pN, our model predicts a relatively narrow window of critical damage stretch thresholds, from 1.01 to 1.30, which agrees well with experimental observations. Our biophysical damage model can help explain the development and progression of axonal damage across the scales and will provide useful guidelines to identify critical damage level thresholds in response to excessive physical forces.

1. Introduction

Brain injury is a major cause of disability and death that is often triggered by an external impact to the head (Hyder et al., 2007; Taylor, 2017). This impact can consist of a single, severe event that immediately leads to traumatic brain injury, or of repeated mild events that gradually result in chronic traumatic encephalopathy. In both cases, the effect of the impact manifests itself at a much smaller scale in the brain: the scale of the axon (Johnson et al., 2013; Smith and Meaney, 2016).

The axon is part of the nerve cell, the neuron, that further consists of a cell body with the cell nucleus and synapses that form connections with other neurons. Figure 1 illustrates a typical axon as a long and slender protrusion from the cell body to provide signaling pathways and transport highways within and away from the brain (Debanne et al., 2011). The axonal cytoskeleton consists of a system of longitudinally aligned microtubules and neurofilaments (Ouyang et al., 2013; Kirkcaldie and Collins, 2016) surrounded by an actin cortex and layers of fatty material, the myelin sheath. Axonal microtubules are composed of heterodimers of α- and β-tubulin that form 13 laterally joined protofilaments, each up to 100μm long (Alberts et al., 2014). These microtubules never run continuously from the cell body to the distal end of the axon. Instead, they form overlapping segments with 10–50 microtubules in any given cross section (Krieg et al., 2017). Microtubules are interconnected by active and passive crosslinking proteins including dynein, kinesin, and tau (Coles and Bradke, 2015). Recent studies suggest that tau protein plays a major role in various types of neurodegeneration that are collectively recognized as tauopathies. Classical examples include Alzheimer's disease, Pick's disease, progressive supranuclear palsy, and chronic traumatic encephalopathy (Woerman et al., 2016; Eisenberg and Sawaya, 2017). A classical hallmark of chronic traumatic encephalopathy is an abnormal increase of tangled tau protein across the brain (Mez et al., 2017). Yet, the precise cause, development, and diagnosis of chronic traumatic encephalopathy are only incompletely understood and remain active fields of research (Asken et al., 2017).

FIGURE 1

Figure 1. The axon is a long and slender protrusion from the neuronal cell body that consists of a system of longitudinally aligned microtubules. Microtubules are composed of heterodimers of α- and β-tubulin, shown in green and blue, that form 13 laterally joined protofilaments, each up to 100 μm long. Axons can extend several centimeters in length and their microtubules never run continuously from the cell body to the distal end. Instead, they form overlapping segments with 10–50 microtubules in any given cross section. Neuronal microtubules are stabilized and cross-linked by tau proteins, shown in red, which bind to microtubules with their three or four binding repeats, shown in yellow. Several components of the axon including neurofilaments, other crosslinking proteins, and cytoskeletal organelles are not displayed.

Physical forces play an important role in the axon under physiological conditions (Suter and Miller, 2011; O'Toole et al., 2015). However, beyond a critical level, forces can trigger axonal degradation and damage (van den Bedem and Kuhl, 2015). Indeed, physical impacts to the head that result in excessive axonal stretch (Ji et al., 2014) may trigger a gradual degradation of the tau-microtubule complex (van den Bedem and Kuhl, 2017). Tau protein is an intrinsically disordered protein with three or four binding repeats that bind to neuronal microtubules and prevents them from depolymerization (Kadavath et al., 2015). Bound tau protein is believed to form an electrostatic zipper with tau protein from neighboring microtubules (Fitzpatrick et al., 2017). As such, tau plays a critical role in stabilizing individual microtubules (Chung et al., 2015) and forming microtubule bundles (Choi et al., 2016). Increasing evidence suggests that axonal damage develops when a physical force is large enough to break the tau-microtubule bonds. An excessive loss of tau crosslinks results in the depolymerization of microtubules (Kadavath et al., 2015), the disintegration of microtubule bundles (Krieg et al., 2017), and the disruption of axonal transport (Tang-Schomer et al., 2012). While experimental testing of the structural integrity of the tau-microtubule complex remains challenging (Li et al., 2015), computational modeling of the axon in response to physical forces can provide useful mechanistic insight into these causal relations (de Rooij and Kuhl, 2018). Although all cytoskeletal elements contribute to the mechanical properties of the axon (Kirkcaldie and Collins, 2016), recent studies suggest that the mechanical stiffness of the axon is most reduced when disrupting axonal microtubules, followed by neurofilaments and microfilaments (Ouyang et al., 2013). Mechanical models of the axon have therefore mainly focused on microtubules (Suter and Miller, 2011), which, because of their hollow circular cross section, provide the largest resistance to bending and tension (Howard, 2001). Models of the axon exist at various levels of complexity ranging from a combination of rheological spring and dashpot elements (O'Toole et al., 2008) via a discrete arrangement of microtubules and crosslinks (Peter and Mofrad, 2012; Ahmadzadeh et al., 2015; Jakobs et al., 2015; Lazarus et al., 2015) to a continuum representation of the axon as a whole (Recho et al., 2016; García-Grajales et al., 2017). Our group has recently proposed a new axonal damage model that integrates the dynamics of microtubule polymerization and depolymerization, the biology of crosslink attachment and detachment, and physics of stretching using a custom-designed finite element model (de Rooij et al., 2016; de Rooij and Kuhl, 2018).

Although brain damage has its mechanistic origin at the axon level, the severity of a head impact is often quantified at the whole brain level by applying experimentally measured linear and rotational accelerations to a human head model (Kuo et al., 2017). In these models, it is essential to accurately capture the brain geometry (Kleiven and von Holst, 2002; Takhounts et al., 2003) and its material properties (de Rooij and Kuhl, 2016; Budday et al., 2017). The most critical input to these models, however, is the critical strain or stretch level beyond which axonal damage occurs (Bain et al., 2004). To better understand the propagation of axonal damage across the scales, we have to connect the axon level to the whole brain level. Toward this goal, we simulate the effect of physical forces across the axon and derive a continuum model for axonal damage as function of the applied stretch and stretch rate. Central to our model is the classical Bell model (Bell, 1978) that characterizes the dynamics of the tau-microtubule complex, from which we infer a damage evolution law that can be easily embedded into finite element models at the whole brain level. Our work provides a systematic strategy to mechanistically correlate crosslink dynamics on the microscopic scale to the evolution of axonal damage on the mesoscopic scale. We anticipate that this work will provide insight into the development of brain damage across the scales and improve current modeling techniques to quantify brain damage for a given physical impact to the brain.

2. Methods

2.1. Axon Model

We model of the axon as a system of straight microtubules that are aligned in the longitudinal direction. Each cross section of the axon has 19 potential microtubules sites arranged in a triangular grid (de Rooij et al., 2016). On average, only half of these potential sites are occupied by a microtubule. As Figure 2 indicates, we assume that all microtubules have the same length and are randomly distributed along the axon. In our finite element model, each microtubule consists of 1,250 one-dimensional truss elements (de Rooij and Kuhl, 2018).

FIGURE 2

Figure 2. Axon model with longitudinally aligned microtubules that are connected by tau protein crosslinks. To account for axonal dynamics, we model each crosslink as a noncovalent slip bond and assign each crosslink a mechanism of Bell model type. We fix the axon at its distal, left end and apply a stretch, λ, and stretch rate, $\stackrel{·}{\lambda }$, to its proximal, right end.

Neighboring microtubules within a cross section are crosslinked by tau protein. At the beginning of the simulation, these crosslinks are randomly distributed across the axon based on a crosslink density parameter. To account for the dynamic behavior of the axon, we have created an extension to the standard finite element method that can either effect individual finite elements or to sets of finite elements (de Rooij et al., 2016). This dynamic behavior represents the molecular mechanisms of particular proteins. Here, we assign a mechanism to each tau protein crosslink. The mechanism describes the dynamic behavior of crosslink detachment and reattachment by modeling the crosslink as a slip bond. The detachment and attachment rates, k(F), are governed by the classical Bell model (Bell, 1978) that characterizes the strength of a chemical bond under a mechanical force F:

$\begin{array}{cc}k(F)=\{\begin{array}{cc}{k}_0\hfill & \text{attach}\\ k_0\exp (F/F_0)& \text{detach}\end{array}\text{ with }{F}_0=\frac{k_{B}T}{\xi }\hfill & (1)\end{array}$

where k₀ is the rate of crosslink attachment and detachment due to random thermal fluctuations. According to the Bell model for slip bonds, the likelihood of detachment increases exponentially with the physical force F applied to the bond. The sensitivity to this mechanical force is described by the characteristic bond strength F₀ = k_BT/ξ, where k_B is the Boltzmann constant, T is temperature, and ξ is the characteristic bond separation distance.

We fix the axon at its proximal end, in our model on the left side, where it connects to the cell body with the nucleus, and apply a physical stretch at the distal end, on the right side, where it connects to other axons. This implies that we apply homogeneous Dirichlet boundary conditions to the microtubules at both ends, zero on the left and non-zero on the right (de Rooij et al., 2016). We constrain all nodes in the model to move along the longitudinal axon direction only. To represent the axonal cytosol, we embed the axon in a viscous fluid with a viscosity of 5 mPa·s (Haak et al., 1976). Table 1 provides an overview of all model parameters.

TABLE 1

Table 1. Parameters of the axon model, the microtubule model, and the crosslink model.

Figure 3 illustrates the flowchart to solve our discrete axon model within our custom-designed finite element algorithm. To include the dynamic behavior of individual proteins, we extended the conventional finite element method by mechanisms (de Rooij et al., 2016). We include a mechanism by introducing a Mechanism object with an Apply() function that contains the full description of the crosslink behavior. We assign this mechanism to each crosslink and execute the Apply() function at the beginning of each iteration step in the Newton-Raphson solver with adaptive time stepping. To apply the slip bond mechanism of the Bell model, the Apply() function has to ensure that each crosslink detaches and attaches at a rate k(F) as in Equation (1) (de Rooij and Kuhl, 2018). At each time step of our simulation, we calculate the probability of detachment or reattachment. We distinguish two cases to calculate the probability of detachment: detachment at a constant force F and at a linearly increasing force F = r_ft.

FIGURE 3

Figure 3. Axon model flowchart based on a conventional Newton-Raphson solver (blue) modularly enhanced by the application of mechanisms (orange). At the beginning of each Newton-Raphson step, we apply the mechanism to each crosslink element. In case the algorithms does not converge, we restore the last equilibrium state of all mechanisms and reduce the time step size.

2.1.1. Crosslink Detachment at a Constant Force F

For a constant force F, we compute the probability of crosslink detachment, p(F, t), at time t, based on the detachment rate k(F) in Equation (1):

$\begin{array}{rc}p(F,t)=k(F)exp(-k(F)t).& (2)\end{array}$

This probability function directly yields the probability of detachment within the next time step, between t₀ and t₀ + Δt, as:

$\begin{array}{cc}P(F,\Delta t)=\frac{{\displaystyle {\int }_{t_0}^{t_0+\Delta t}p}(F,\overline{t})d\overline{t}}{{\displaystyle {\int }_{t_0}^{\infty }p}(F,\overline{t})d\overline{t}}=1-\exp (-k(F)\Delta t).& (3)\end{array}$

To obtain the probability of attachment of a crosslink, we simply substitute F = 0 into Equation (3).

2.1.2. Crosslink Detachment at a Linearly Increasing Force F = r_ft

For a linearly increasing force F, the probability of detachment becomes (de Rooij and Kuhl, 2018):

$\begin{array}{rc}p(F,r_f)=\frac{k(F)}{r_f}exp\left(-\frac{F_0}{r_f}\left[k(F)-k_0\right]\right).& (4)\end{array}$

The probability of detachment within the next time step, between t₀ and t₀+Δt, then follows as:

$\begin{array}{cc}P(F,\Delta F)=\frac{{\displaystyle {\int }_{F_0}^{F_0+\Delta F}p}(F)dF}{{\displaystyle {\int }_{F_0}^{\infty }p}(F)dF},& (5)\end{array}$

where ΔF = r_fΔt. In our simulations, we compute the individual loading rate r_f for each crosslink as the crosslink force divided by the time the crosslink has been attached.

2.2. Discrete Axonal Damage Model

The major objective of our axon model is to interpret axonal damage as an emergent feature of the dynamic attachment and detachment of crosslinks. Motivated by the common definition of damage in continuum damage mechanics (Kachanov, 1986), we define axonal damage as the relative loss of axonal stiffness due to excessive detachment of crosslinks as the result of a physical force. Indeed, Equation (1) shows that a physical force F increases the detachment rate of crosslinks and, thereby, results in net reduction of attached crosslinks. To quantify the amount of damage, we first need to characterize the baseline, undamaged, mechanical response of our model axon. This baseline response is nonlinear, viscoelastic, and, therefore, rate dependent (de Rooij and Kuhl, 2018). We obtain the baseline response by performing simulations at an infinite characteristic bond strength F₀ → ∞, see Equation (1). This implies that the detachment rate, k(F; F₀ → ∞) = k₀, is the same as the attachment rate, and the total number of attached crosslinks will, on average, remain constant. To account for viscous effects in the baseline response of the axon, we perform these baseline simulations for a range of loading rates r_f.

We define axonal damage for a finite characteristic bond strength F₀ as the relative stiffness degradation with respect to the undamaged axon. Consistent with continuum damage mechanics (Lemaitre, 1992), we use the scalar-valued damage parameter d to quantify axonal damage. Damage ranges from d = 0 for a completely intact axon to d = 1 for a fully damaged axon and relates the reduced stiffness E to the initial undamaged stiffness E₀ as:

$\begin{array}{rc}E=\left[1-d\right]E_0,& (6)\end{array}$

where both E and E₀ are the corresponding secant stiffnesses. Figure 4 (left) shows characteristic force-stretch curves of the axon for varying characteristic bond strengths F₀. From these force-stretch curves, we compute axonal damage d(λ) for a given stretch λ as:

$\begin{array}{rc}d(\lambda )=1-E(\lambda )/E_0(\lambda ).& (7)\end{array}$

FIGURE 4

Figure 4. Characteristic force vs. stretch and damage vs. stretch curves for varying characteristic bond strengths, F₀. An infinite bond strength F₀ → ∞ defines the baseline, undamaged stiffness E₀(λ). Finite bond strengths F₀ trigger a net increase in crosslink detachment resulting in an increase in damage d and a reduced axonal stiffness E(λ). Increasing the bond strength F₀ decreases the not crosslink detachment, which increases the stretch λ and force F required to trigger axonal damage. We define axonal damage, d(λ) = 1−E(λ)/E₀(λ), as the ratio of the reduced and undamaged secant stiffnesses E(λ) and E₀(λ).

Figure 4 (right) shows characteristic damage-stretch curves of the axon for varying characteristic bond strengths F₀.

2.3. Homogenization

To bridge the scales, we postulate a specific functional form for the evolution of damage and assume that its parameters emerge naturally from the dynamic behavior of the crosslinking tau proteins. In other words, we seek an evolution equation that provides an analytical expression for axonal damage, $d(\lambda ,\stackrel{·}{\lambda })$, as function of the axonal stretch λ and stretch rate $\stackrel{·}{\lambda }$. We characterize the accumulation of damage through the logistic function (Verhulst, 1845),

$\begin{array}{cc}d(\lambda ,\dot{\lambda })=\widehat{\mathscr{H}}(\lambda -{\lambda }_{50}(\dot{\lambda };F_0);\alpha (\dot{\lambda };F_0)).& (8)\end{array}$

The C^∞-smooth Heaviside function $\widehat{{H}}(x;\alpha )=e^{\alpha x}/\left[1+e^{\alpha x}\right]$ represents an S-shaped sigmoid curve, λ₅₀ is the half damage stretch at the midpoint of the S-shaped curve, at which d(λ = λ₅₀) = 0.5, and α is proportional to the slope at this midpoint. Widely used in population dynamics, the logistic function implies that the initial stage of damage is approximately exponential, it then begins to saturate at the half damage stretch λ₅₀, and gradually converges to the fully damaged state, d = 1, as the stretch increases. We assume that both λ₅₀ and α are functions of the characteristic bond strength F₀ and vary with the applied stretch rate $\stackrel{·}{\lambda }$. Importantly, in our model, this rate dependence emerges naturally from the underlying crosslink dynamics.

2.3.1. Ansatz for the Half Damage Stretch ${\lambda }_{50}$

For the half damage stretch ${\lambda }_{50}(\stackrel{·}{\lambda };F_0)$, we use Equation (4) to compute the expected crosslink force, $\widehat{F}$, at which a crosslink detaches:

$\begin{array}{cc}\widehat{F}={\displaystyle {\int }_0^{\infty }\tilde{F}}p(\tilde{F};r_f)d\tilde{F}=F_0\exp \left(\frac{k_0{F}_0}{r_{\text{}f}}\right)\Gamma \left(0,\frac{k_0{F}_0}{r_{\text{}f}}\right)& (9)\end{array}$

where $\Gamma (a,b)={\int }_b^{\infty }{e}^{-x}{x}^{a-1}dx$ is the upper incomplete gamma function. When the crosslinks are attached to microtubules, we assume that they behave linearly elastic and we expect the transition stretch to be proportional to the detachment force, $\left[{\lambda }_{50}-1\right]\propto \widehat{F}$. With $k_0{F}_0/r_f\propto 1/\stackrel{·}{\lambda }$, we therefore propose:

$\begin{array}{rc}{\lambda }_{50}(\stackrel{·}{\lambda };F_0)=1+a_{\lambda }exp\left(\frac{b_{\lambda }}{\stackrel{·}{\lambda }}\right)\Gamma \left(0,\frac{b_{\lambda }}{\stackrel{·}{\lambda }}\right),& (10)\end{array}$

where a_λ(F₀) and b_λ(F₀) are parameters that depend on the characteristic bond strength F₀ and will emerge naturally from the tau crosslink dynamics within the axon.

2.3.2. Ansatz for the Damage Slope $\alpha$

For the damage slope $\alpha (\stackrel{·}{\lambda };F_0)$, we follow a similar approach. Since we interpret damage as the net loss of crosslinks, axonal damage is a function of the fraction of attached crosslinks ${\widehat{n}}_{\text{att}}$:

$\begin{array}{rc}d=1-{\widehat{n}}_{\text{att}}\text{ with}{\widehat{n}}_{\text{att}}=\frac{{\widehat{t}}_{\text{att}}}{\frac{1}{2}\left[{\widehat{t}}_{\text{att}}+{\widehat{t}}_{\text{det}}\right]}& (11)\end{array}$

where ${\widehat{t}}_{\text{att}}$ and ${\widehat{t}}_{\text{det}}$ are the expected duration that a crosslink is attached and detached respectively. This implies that, in the limit of homeostasis between attachment and detachment, ${\widehat{t}}_{\text{att}}={\widehat{t}}_{\text{det}}$ and ${\widehat{n}}_{\text{att}}=1$ and d = 0, whereas in the limit of an excessive detachment, ${\widehat{t}}_{\text{att}}<<{\widehat{t}}_{\text{det}}$ and ${\widehat{n}}_{\text{att}}=0$ and d = 1. From the definition of the S-curve, we know that α is proportional to the slope of the damage curve at d = 0.5:

$\begin{array}{cc}\alpha \propto {\frac{\text{d}d}{\text{d}\lambda }|}_{\text{d}=0.5}={\frac{\text{d}d}{\text{d}{\widehat{t}}_{\text{att}}}·\frac{\text{d}{\widehat{t}}_{\text{att}}}{\text{d}\dot{\lambda }}·\frac{\text{d}\dot{\lambda }}{\text{d}\lambda }|}_{\text{d}=0.5}\propto \frac{\text{d}{\widehat{t}}_{\text{att}}}{\text{d}\dot{\lambda }}.& (12)\end{array}$

By combining the definition of the attachment time, ${\widehat{t}}_{\text{att}}=\widehat{F}/r_f$, with Equations (9) and (12), we propose:

$\begin{array}{rc}\alpha (\stackrel{·}{\lambda };F_0)=\frac{a_{\alpha }}{{(\stackrel{·}{\lambda })}^3}\left[\stackrel{·}{\lambda }-\left[b_{\alpha }+\stackrel{·}{\lambda }\right]exp\left(\frac{b_{\alpha }}{\stackrel{·}{\lambda }}\right)\Gamma \left(0,\frac{b_{\alpha }}{\stackrel{·}{\lambda }}\right)\right],& (13)\end{array}$

where, again, a_α(F₀) and b_α(F₀) are parameters that depend on the characteristic bond strength F₀ and will emerge naturally from the tau crosslink dynamics within the axon.

2.4. Continuum Axonal Damage Model

To embed the homogenized equations into a continuum axonal damage model, we introduce the deformation φ(X, t) along the axis of the axon and define the axon level stretch λ and stretch rate $\stackrel{·}{\lambda }$,

$\begin{array}{rc}\lambda =\frac{\partial \phi (X,t)}{\partial X}\text{\hspace{1em}and\hspace{1em}}\stackrel{·}{\lambda }=\frac{\mathrm{\text{d}}\lambda (X,t)}{\mathrm{\text{d}}t}.& (14)\end{array}$

We then introduce the free energy density of the damaged axon Ψ as the damage weighted stored energy of the undamaged, elastic axon Ψ₀,

$\begin{array}{cc}\Psi (\lambda ,\dot{\lambda })=[1-d]{\Psi }_0(\lambda )\text{\hspace{1em}with\hspace{1em}}d=d(\lambda ,\dot{\lambda }),& (15)\end{array}$

and assume that the evolution of damage is driven by both stretch and stretch rate, $d(\lambda ,\stackrel{·}{\lambda })$, while the elastic energy is a function of the stretch alone Ψ₀(λ). Motivated by standard arguments of thermodynamics, we introduce the Cauchy stress σ = Pλ and the Piola stress P as thermodynamically conjugate quantity to the stretch λ, and interpret the Piola stretch P as damage weighted elastic axonal stress P₀,

$\begin{array}{rc}P=\frac{\partial \Psi }{\partial \lambda }=\left[1-d\right]P_0\text{\hspace{1em}with\hspace{1em}}{P}_0=\frac{\partial {\Psi }_0}{\partial \lambda }.& (16)\end{array}$

To keep track of the maximum amount of stretch the axon has experienced throughout its history, it is common practice to introduce an internal variable,

$\begin{array}{rc}{\lambda }^{*}={\displaystyle \underset{0\le t\le \tau }{max}}\left\{\lambda (t)\right\},& (17)\end{array}$

which drives the evolution of damage,

$\begin{array}{rc}d=\frac{exp(\alpha \left[{\lambda }^{*}-{\lambda }_{50}\right])}{1+exp(\alpha \left[{\lambda }^{*}-{\lambda }_{50}\right])}.& (18)\end{array}$

The stretch rate dependent half damage stretch,

$\begin{array}{rc}{\lambda }_{50}(\stackrel{·}{\lambda };F_0)=1+a_{\lambda }exp\left(\frac{b_{\lambda }}{\stackrel{·}{\lambda }}\right)\Gamma \left(0,\frac{b_{\lambda }}{\stackrel{·}{\lambda }}\right),& (19)\end{array}$

and the stretch rate dependent damage slope,

$\begin{array}{rc}\alpha (\stackrel{·}{\lambda };F_0)=\frac{a_{\alpha }}{{(\stackrel{·}{\lambda })}^3}\left[\stackrel{·}{\lambda }-\left[b_{\alpha }+\stackrel{·}{\lambda }\right]exp\left(\frac{b_{\alpha }}{\stackrel{·}{\lambda }}\right)\Gamma \left(0,\frac{b_{\alpha }}{\stackrel{·}{\lambda }}\right)\right],& (20)\end{array}$

follow from the homogenization in section 2.3 and vary with the bond strength F₀ of the individual crosslinks in the axon. Last, to solve the continuum equations of axonal damage within a finite element setting, we derive the tangent modulus,

$\begin{array}{rc}\mathrm{\text{A}}=\frac{\mathrm{\text{d}}P}{\mathrm{\text{d}}\lambda }=\left[1-d\right]{\mathrm{\text{A}}}_0-\frac{\mathrm{\text{d}}d}{\mathrm{\text{d}}\lambda }{P}_0\text{\hspace{1em}with\hspace{1em}}{\mathrm{\text{A}}}_0=\frac{\mathrm{\text{d}}{P}_0}{\mathrm{\text{d}}\lambda }.& (21)\end{array}$

For our specific damage model with

$\begin{array}{ll}\hfill & \frac{\text{d}d}{\text{d}\lambda }=\frac{\text{d}d}{\text{d}{\lambda }^{*}}\frac{\text{d}{\lambda }^{*}}{\text{d}\lambda }\text{\hspace{1em}with\hspace{1em}}\frac{\text{d}d}{\text{d}{\lambda }^{*}}=\alpha [1-d]d\hfill \\ \hfill & \text{and\hspace{1em}}\frac{\text{d}{\lambda }^{*}}{\text{d}\lambda }=\{\begin{array}{l}\mathrm{1...}\text{loading}\hfill \\ \mathrm{0...}\text{unloading}\hfill \end{array}\hfill & (22)\end{array}$

we obtain the following simple structure of the tangent modulus,

$\begin{array}{rc}\text{A}=\left[1-d\right]\left[{\text{A}}_0-\alpha d{P}_0\right].& (23)\end{array}$

For example, for an elastic material of Mooney Rivlin type, with ${\Psi }_0=c_1\left[{\lambda }^2+2/\lambda -3\right]+c_2\left[2\lambda +1/{\lambda }^2-3\right]$, the Cauchy stress becomes $\sigma =\left[1-d\right]2\left[c_1+c_2/\lambda \right]\left[{\lambda }^2-1/\lambda \right]$, the elastic tangent modulus is ${\mathrm{\text{A}}}_0=2c_1\left[1+2/{\lambda }^3\right]+6c_2/{\lambda }^4$, and the elastic Piola stress is $P_0=2\left[c_1+c_2/\lambda \right]\left[\lambda -1/{\lambda }^2\right]$, where $c_1+c_2=\frac{1}{2}\mu$ are the two constitutive parameters of the Mooney Rivlin model and μ is the overall shear modulus of the axon (Goriely et al., 2015a).

3. Results

3.1. Axon Model

Figure 5 shows the result of a single simulation in which we applied a stretch of λ = 1.2 at a stretch rate $\stackrel{·}{\lambda }=10$/s assuming a bond strength of F₀ = 5pN.

FIGURE 5

Figure 5. Characteristic output of a single simulation. The kymograph, (Top left), tracks the position of all microtubules throughout the simulation. It reveals a complete separation between the proximal and distal ends of the axon. This separation manifests itself in a loss of the required external force and an increase in axonal damage (Bottom left). Four snapshots of the axon show the crosslinks color coded by crosslink stretch (Right). These snapshots show an initial increase in crosslink stretch, followed by a decrease after the proximal and distal sets of microtubule bundles have separated.

The computational kymograph in the top left traces the location of all microtubules in the axon throughout the entire duration of the simulation. The snapshots on the top right show the axon at four different time points with the tau crosslinks color coded by the crosslink stretch. Toward the end of the simulation, the stretch in the remaining crosslinks decreases. This decrease is accompanied by a clear separation between the proximal and distal ends of the axon, as we can see from the kymograph. This separation is indicative of axonal damage; it reduces the overlap of microtubules and, thereby, the number of connecting crosslinks in the damaged region. Indeed, the force vs. time curve in the bottom left shows an initial increase in force followed by a rapid decrease approximately 5 ms into the simulation. The axonal damage, computed according to section 2.2, increases from d = 0 at the beginning of the simulation to d = 1 when the two sets of microtubule bundles are fully disconnected.

3.2. Discrete Axonal Damage

To probe the sensitivity of axonal damage with respect to stretch and stretch rates, we performed several sets of simulations for a range of stretches, λ∈[1, 1.2], and stretch rates, $\stackrel{·}{\lambda }\in \left[0.1,10\right]$ /s. The only input parameter that is not well defined in the literature is the characteristic bond strength of the tau crosslinks. We therefore also probed a range of crosslink bond strengths, F₀ = [1, 5, 10, 50, 100] pN.

Figure 6 (top) shows the damage vs. stretch curves for n = 1,000 simulations at all five crosslink bond strengths, F₀ = [1, 5, 10, 50, 100] pN, color coded by the stretch rate $\stackrel{·}{\lambda }\in \left[0.1,10\right]$ /s. Consistent with the Bell model, for smaller values of F₀, axonal damage develops earlier, at lower axonal stretch levels, than for larger values of F₀. Indeed, Equation (1) shows that at lower crosslink bond strengths F₀, the detachment rates for a given applied crosslink force F are higher, which manifests itself in an increased axonal damage. Our axon model also predicts that crosslink detachment is more likely to happen at higher crosslink stretches for high loading rates, which follows directly from Equation (9). Figure 6 captures this prediction as axonal damage develops at higher values for axonal stretch for high loading rates. Interestingly, this trend is reversed when considering damage vs. loading time.

FIGURE 6

Figure 6. Axonal damage vs. axonal stretch color coded by the stretch rate for five bond strengths. Discrete axonal damage simulations with n = 10 simulations for each stretch rate (Top); homogenization using the average response for each stretch rate to identify the parameters for S-shaped curve of the homogenized damage model (Middle); and continuum axonal damage simulation (Bottom). The graphs demonstrate an important feature of the Bell model: at higher stretch rates, axonal damage occurs at higher axonal stretch.

Figure 7 shows the same simulations as Figure 6, but now as damage plotted vs. loading time. It is clearly visible that higher applied stretch rates lead to earlier development of axonal damage. Thus, increased loading rates triggers earlier development of axonal damage, but at a higher axonal stretch, all consistent with the Bell model.

FIGURE 7

Figure 7. Axonal damage vs. time color coded by the stretch rate for five bond strengths. Discrete axonal damage simulations with n = 10 simulations for each stretch rate (Top); homogenization using the average response for each stretch rate to identify the parameters for S-shaped curve of the homogenized damage model (Middle); and continuum axonal damage simulation (Bottom). The graphs represent the same data as in Figure 6, and demonstrate an important feature of the Bell model: at higher stretch rates, axonal damage occurs at higher axonal stretch, but earlier in time.

3.3. Homogenization

To homogenize the results or our discrete axon model simulation toward an overall constitutive damage model for the axon, we calibrate our damage model $d(\lambda ,\stackrel{·}{\lambda })$ in Equations (8), (19), and (20), using our discrete axon level simulations. For each characteristic crosslink force and applied stretch rate, we compute the mean damage vs. stretch curve from our simulations. For each mean curve, we calibrate the half damage stretch λ₅₀ and the damage slope α according to Equation (8). Figures 6, 7 (middle) illustrate the mean damage curves for each applied stretch rate, together with the best analytical fit as a dashed, black line. These figures show a good qualitative agreement between the discrete axon model and its homogenized response, which supports our initial selection for the damage evolution law, $d(\lambda ,\stackrel{·}{\lambda })$, in Equation (8).

For each applied stretch rate, $\stackrel{·}{\lambda }$, the homogenization yields one value for the half damage stretch λ₅₀ and for the damage slope α, assuming a fixed characteristic bond strength, F₀. In the next step, we use these values together with Equations (19) and (20) to obtain a stretch-rate dependent half damage stretch ${\lambda }_{50}(\stackrel{·}{\lambda })$ and damage slope $\alpha (\stackrel{·}{\lambda })$. From the best fits, we obtain discrete values for the parameters a_λ and b_λ in Equation (19) and a_α and b_α in Equation (20). Naturally, these four parameters will be different for different characteristic bond strengths F₀.

Figure 8 shows the numerical data points and the analytical fits for the half damage stretch ${\lambda }_{50}(\stackrel{·}{\lambda })$ and for the damage slope $\alpha (\stackrel{·}{\lambda })$ for the range of F₀ = [1, 5, 10, 50, 100] pN. Qualitatively, Figure 8 confirms that our expressions for ${\lambda }_{50}(\stackrel{·}{\lambda };F_0)$ and $\alpha (\stackrel{·}{\lambda };F_0)$ in Equations (19) and (20), respectively, accurately represent the simulation data. In addition, Figure 8 confirms that λ₅₀ increases with increasing stretch rate, $\stackrel{·}{\lambda }$, and with increasing crosslink bond strength, F₀. Both trends are consistent with our crosslink model and with Figure 6. In contrast, the damage slope parameter α decreases with increasing stretch rate, $\stackrel{·}{\lambda }$, and with increasing crosslink bond strength, F₀.

FIGURE 8

Figure 8. Homogenization of half damage stretch λ₅₀ and damage slope α from average response of n = 10 discrete axonal damage simulations for as function of the stretch rate, $\stackrel{·}{\lambda }$, for five bond strengths, F₀. Blue dots represents λ₅₀ and α values from the average damage vs. stretch relations in Figure 6; solid black lines represent homogenization using Equations (19) and (20) for five bond strengths, F₀.

Table 2 summarizes the homogenized parameters a_λ, b_λ, a_α, and b_α, for bond strengths within the range F₀ = [1, 5, 10, 50, 100] pN. With increasing bond strength F₀, a_λ, b_λ, and b_α increase, while a_α decreases. To interpolate between the five bond strengths, we suggest the following rationale: Motivated by Equations (19) and (20), we expect that a_λ → 0 and a_α → ∞ as the characteristic bond strength decreases toward zero, F₀ → 0. This suggests power law relations for the parameters a of the form $a_{\lambda }=1.936·10^{-3}{(F_0)}^{0.835}$ and $a_{\alpha }=-0.675{(F_0)}^{-0.166}$. Motivated by phenomenological considerations, for the parameters b, we suggest a linear dependence on ln(F₀) of the form $b_{\lambda }=4.64·10^{-6}ln(F_0)+1.61·10^{-6}$ and $b_{\alpha }=1.55·10^{-2}ln(F_0)+4.55·10^{-2}$. This completes our damage model that is fully determined as function of stretch, λ, and stretch rate $\stackrel{·}{\lambda }$, parameterized by the characteristic bond strength, F₀.

TABLE 2

Table 2. Homogenized axon parameters a_λ and b_λ to calculate the half damage stretch λ₅₀ according to Equation (19) and a_α and b_α to calculate the damage slope α according to Equation (20) for a range of characteristic bond strengths F₀.

3.4. Continuum Axonal Damage

Once homogenized and calibrated, we can use the axonal damage model and embed it into a continuum damage simulation using the governing equations from section 2.4. We can, for example, embed these equations into a nonlinear finite element analysis and project the homogenized axonal response along the axonal direction, in a one-, two-, or fully three-dimensional brain model. These continuum damage simulations are fully determined by our evolution equations for the damage variables $d(\lambda ,\stackrel{·}{\lambda })$, which emerge naturally from the axon-level crosslink dynamics.

Figures 6, 7 (bottom) summarize the resulting damage vs. stretch and damage vs. time contours color coded by the stretch rate for five bond strengths. Both graphs highlight an important feature of the Bell model: at higher stretch rates, axonal damage occurs at higher axonal stretch, but earlier in time. Comparing the continuum axonal damage model (bottom) to the discrete axonal damage model (top) and its homogenization (middle) confirms that our transient damage model at the continuum level captures the same failure characteristics as the discrete axonal damage model based on crosslink detachment and reattachment dynamics.

Figure 9 illustrates the emergent axonal damage vs. stretch λ[1.0, 1.5] (left) and time t∈[0.0, 1.0] s (right) at varying stretch rates $\stackrel{·}{\lambda }\in \left[0.1,10\right]$/s at a constant bond strength of F₀ = 100 pN. The white circles represent experimentally characterized damaged, d = 1, and undamaged, d = 0, nervous tissue samples that had been exposed to different strain levels (Bain and Meaney, 2000). These experiments clearly report a transition from a low likelihood of damage at low stretch rates to a high likelihood of damage at high stretch rates and characterize the critical stretch levels at which damage emerges. The black circles define the conservative damage threshold at d = 0.05 for 14% strain and a stretch of 1.13, the liberal damage threshold at d = 0.90 for 34% strain and a stretch of 1.30, and the optimal damage threshold at d = 0.25 for 21% strain and a stretch of 1.19, indicated through the black dashed dlines (Bain and Meaney, 2000). The thick black line highlights the best fit. Its half damage stretch is λ₅₀ = 1.22, which implies that at a stretch of 1.22 corresponding to a 25% strain, half of the samples were damaged. The graphs in Figure 9 demonstrate an emergent feature of our transient crosslink model: For increasing stretch rates, damage develops at a larger stretch (left) but earlier in time (right).

FIGURE 9

Figure 9. Continuum axonal damage. Axonal damage vs. axonal stretch (Left) and time (Right) at varying stretch rates and constant bond strength. White circles represent experimentally characterized damaged, d = 1, and undamaged, d = 0, nervous tissue samples at different strain levels (Bain and Meaney, 2000). Black circles define the conservative damage threshold at 14% strain, the liberal damage threshold at 34% strain, and the optimal damage threshold at 21% strain (Bain and Meaney, 2000). The graphs demonstrate an emergent feature of our transient crosslink model: For increasing stretch rates, damage develops at a larger stretch (Left) but earlier in time (Right).

Figure 10 illustrates the emergent axonal force vs. stretch λ∈[1.0, 1.5] (left) and time t∈[0.0, 1.0] s (right) at varying stretch rates $\stackrel{·}{\lambda }\in \left[0.1,10\right]$/s at a constant bond strength of F₀ = 100 pN. The dark red area marks the elastic, undamaged regime, here for the example of a Mooney Rivlin model according to section 2.4, with a parameter ratio of c₁:c₂ = 3:1. All other colors highlight the effect of damage with a gradually increasing force that eventually reaches a peak and decreases as a result of axonal softening. The continuum level force vs. time behavior (right) agrees well with the axon level force vs. time behavior in Figure 5 (bottom left). The graphs in Figure 10 demonstrate an emergent feature of our transient crosslink model: For increasing stretch rates, the peak axonal force increases (left) but is reached earlier in time (right).

FIGURE 10

Figure 10. Continuum axonal damage. Axonal force vs. axonal stretch (Left) and time (Right) at varying stretch rates and constant bond strength. The dark red area marks the elastic, undamaged regime. The graphs demonstrate an emergent feature of our transient crosslink model: For increasing stretch rates, the peak axonal force increases (Left) but is reached earlier in time (Right).

4. Discussion

Brain damage can be caused by a wide variety of physical impacts, ranging from a single and strong blow to the head to several mild but repeated concussive events. Independent of the type of impact, brain damage typically originates at the level of the axon: Diffuse axonal injury may develop instantaneously upon a single severe impact, whereas chronic traumatic encephalopathy develops gradually in response to repeated mild impacts to the head. Yet, the precise mechanisms how physical impacts to the head triggers pathologies at the axon level remain incompletely understood (Goriely et al., 2015b). With current technologies, we cannot reliably measure the direct effects of physical forces to the head. However, mechanical modeling can help to indirectly assess the effects of physical impact and correlate external loading to critical damage thresholds on the axonal level (Greenwald et al., 2008; Giordano et al., 2017; Kuo et al., 2017). Toward this goal we propose a mechanistic biophysical model that interprets axonal damage as an emergent property of crosslink dynamics, physical stretches, and stretch rates.

4.1. Axonal Damage Is a Result of Excessive Crosslink Detachment

We model the axon as a parallel arrangement of longitudinally aligned microtubules that are crosslinked by tau protein. These crosslinks can break and form according to the Bell model for chemical bond breaking (Bell, 1978) under external physical forces (Evans and Ritchie, 1997). The Bell model is characterized by a characteristic bond strength, F₀, which is the only unknown variable in our axon model. All other parameters in the model have been reported in the literature as summarized in Table 1. To investigate the evolution of damage inside the axon, we apply a displacement controlled external stretch at different stretch rates. Damage emerges naturally as a result of the stretch-induced forces acting on the crosslinks, which, according to the Bell model, triggers a net increase of crosslink detachment. From a physics perspective, we define axonal damage as the loss in axon stiffness (Kachanov, 1986; Lemaitre, 1992) triggered by a gradual loss of crosslinks (Ahmadzadeh et al., 2014) that promotes microtubule depolymerization (Kadavath et al., 2015) and destabilizes the axonal cytoskeleton and (Krieg et al., 2017).

4.2. Damage Accumulates at the Location of Weakest Connectivity

A representative simulation of axonal damage is characterized by an applied stretch, λ, at a given stretch rate, $\stackrel{·}{\lambda }$. The main output of our simulation is the overall force-stretch behavior of the axon, see Figure 5. We use this axonal force-stretch response to derive the axonal damage-stretch response compared to the undamaged baseline case, see section 2.2. Figure 5 shows that the axonal force increases initially, peaks, and then quickly drops down to zero. This drop is a defining feature of axonal damage and indicates that the axon has lost all its mechanical resistance to loading. The kymograph in Figure 5 illustrates that this rapid loss in stiffness is associated with a primary axotomy, the development of two disconnected sets of microtubule bundles, one at the proximal and one at the distal end, which ultimately defines axonal failure. The exact location of the axotomy is stochastic due to the probabilistic nature of the slip bond model and the underlying axon geometry: Once a weak cross section randomly develops along the axon, each remaining crosslinks in this cross section has to carry more mechanical load; this increases its probability of detachment, which increases the probability that the cross section becomes even weaker and eventually fails completely.

4.3. At Higher Stretch Rates, Axons Can Sustain Higher Stretches

We perform numerical simulations for a range of stretches λ∈[1, 1.2] and stretch rates $\stackrel{·}{\lambda }\in$ [0.1,10] /s, and systematically vary the characteristic bond strength over two orders of magnitude F₀∈ [1,100] pN. Motivated by the randomness in the precise axonal geometry and in the time of detachment and attachment of individual crosslinks, we perform n = 10 simulations for each set of input parameters and use the average result of those simulations for further analysis. Figures 6, 7 (top) show the axonal damage vs. stretch and time for a range of stretch rates. The results in these figures are consistent with an important emergent feature of our transient crosslink model: at higher applied stretch rates, axonal damage develops at a higher axonal stretch, but earlier in time (de Rooij and Kuhl, 2018).

4.4. The Homogenized Axonal Damage Behavior Displays an S-Shaped Response

To derive a constitutive model for axonal damage that we can embed into whole brain damage simulations (Goriely et al., 2015a), we homogenize the discrete axonal response. For each combination of stretch, stretch rate, and characteristic bond strength, we perform n = 10 discrete axonal damage simulations and average their damage-stretch response. We homogenize the discrete model by fitting an S-shaped damage curve through the average damage-stretch response. The S-shaped damage model is similar to damage evolution laws for soft materials that exponentially approach complete damage at d = 1 (Beatty and Krishnaswamy, 2000; Weisbecker et al., 2012). However, several damage evolution laws assume that damage only starts beyond a certain stretch threshold (Peerlings et al., 2001). Here, we choose the smooth S-shaped curve because our crosslink model is also smooth (Bell, 1978) and because it nicely captures the homogenized axonal response in Figures 6, 7 (middle).

4.5. The Damage-Stretch Behavior of the Axon Is Defined by Two Parameters

Our S-shaped damage curve is defined by two parameters with a clear physical interpretation: the half damage stretch λ₅₀ that defines the stretch at which the axon is half damaged and the damage slope α associated with the slope at this half damage stretch. Note that the half damage stretch λ₅₀ is conceptually similar to the damage stretch threshold that has been proposed in literature (Marini et al., 2012; Li, 2016). We assume that both the half damage stretch and the damage slope depend on stretch rate and characteristic crosslink force, ${\lambda }_{50}(\stackrel{·}{\lambda };F_0)$ and $\alpha (\stackrel{·}{\lambda };F_0)$. We derive the qualitative dependence on the stretch rate $\stackrel{·}{\lambda }$ from the Bell model in Equations (19) and (20), which we compare to our discrete axon simulations in Figure 8. Figure 8 reveals that the half damage stretch, λ₅₀, increases with increasing stretch rate and with increasing bond strength. The increase of λ₅₀ with increasing stretch rate $\stackrel{·}{\lambda }$ is consistent with the Bell model that assumes that, at higher stretch rates, the axon can sustain higher stretches prior to damaging. The increase of λ₅₀ with increasing bond strength F₀ is also consistent with the Bell model in Equation (1), since higher bond strengths F₀ require higher forces F to generate the same detachment rate. Figure 8 shows that the slope, α, decreases with increasing stretch rate and increasing bond strength. This is consistent with the Bell model in Equation (20) and with the smoothness of the Bell model. Table 2 summarizes the homogenized axon parameters a_λ and b_λ to calculate the half damage stretch λ₅₀ and the homogenized axon parameters a_α and b_α to calculate the damage slope α for a range of bond strengths F₀. Table 2 suggests that, with increasing bond strength F₀, a_λ, b_λ, and b_α increase while a_α decreases.

4.6. Our Damage Model Agrees Well With Experimental Findings

Section 2.4 and Figures 6, 7 (bottom) summarize our constitutive model for axonal damage. Our damage model is completely determined as a function of stretch λ and stretch rate $\stackrel{·}{\lambda }$, parameterized in the bond strength F₀. Although the characteristic bond strength F₀ is currently unknown, our selected range F₀∈ [1,100] pN lies well within the range of physiological force levels observed at the protein level. For example, dynein protein generates a force of about 1 pN (Mallik et al., 2004), microtubule assembly may generate pushing forces up to 3–4 pN (Dogterom and Yurke, 1997), and the growth cone of the axon generates a total pulling force of about 1–20 nN (Rajagopalan et al., 2010; Hyland et al., 2014). Strikingly, Figure 8 suggests that the values for the half damage stretch λ₅₀ all lie within λ₅₀∈ [1.01,1.3] for our entire range of stretch rates $\stackrel{·}{\lambda }\in$ [0.1,10] /s and bond strengths F₀∈ [1,100] pN. By definition, λ₅₀ is the axonal stretch at which 50% of the axon is damaged. This suggests that we can use the λ₅₀ value as a surrogate measure for the axonal damage level threshold. The range of λ₅₀∈ [1.01,1.30] agrees well with reported critical stretch values for axonal injury: critical values between 1.05 and 1.10 have been found based on animal and physical studies (Margulies and Thibault, 1992), a critical stretch of 1.05 caused mild injury in cortical axons in culture (Yuen et al., 2009), axonal stretch between 1.09 and 1.16 led to axonal injury in rats (Singh et al., 2016), and critical stretches between 1.14 and 1.34 were found at stretch rates between 30 and 60 /s in white matter brain tissue (Bain and Meaney, 2000). Notably, these reported critical damage stretch thresholds are all based on a single severe loading of the axon and mimic the event of traumatic brain injury. To date, there is no reliable data on critical damage stretch thresholds for multiple mild loading of the axon that would mimic the event of chronic traumatic encephalopahy (Asken et al., 2017).

4.7. Our Damage Model Integrates Well Into Finite Element Algorithms

Once calibrated and validated, we can embed our constitutive model for axonal damage in a continuum mechanics model for whole brain simulations and superpose it to the isotropic behavior of the tissue (Cloots et al., 2013; Mao et al., 2013; Goriely et al., 2015b; Giordano et al., 2017; Weickenmeier et al., 2017). Figures 9, 10 show that the continuum implementation of the damage model, e.g., within a nonlinear finite element setting, correctly reproduces the axon level features of damage. Notably, a well-known problem with continuum damage models is that, in the softening regime, the governing equations become ill posed and the numerical solutions become mesh dependent. These issues can be resolved with appropriate regularization techniques (de Borst et al., 1993; Kuhl et al., 2000). A natural regularization technique is to account for the rate dependent nature of damage (Geers et al., 1994). Although we do not explicitly investigate regularization here, the inherent rate dependence of our axonal damage model potentially regularizes the simulation at no additional cost (Pereira et al., 2017).

4.8. Our Axon Model Has a Few Limitations

We have proposed a consistent strategy to relate microscale protein behavior to axonal damage and to develop a constitutive damage model that can be used at the continuum, whole brain level. However, we do recognize several limitations to our model that we plan to address in future work: First, our model assumes that axonal damage is solely caused by the disruption in tau protein crosslinks. Although this disruption is consistent with, e.g., the diagnoses of chronic traumatic encephalopathy through an abnormal increase of unbound tau, other mechanisms such as microtubule rupture may contribute to damage of the axon (Tang-Schomer et al., 2010; van den Bedem and Kuhl, 2015). In addition, we simplified the axonal cytoskeleton as a composition of microtubules and tau proteins, while we recognize that the axon contains additional cytoskeletal elements and organelles that could be structurally relevant, such as neurofilaments, microfilaments, dynein, myosin, and the actin cortex (Ouyang et al., 2013; Kirkcaldie and Collins, 2016; Tofangchi et al., 2016). For example, neurofilaments and microfilaments contribute to the axon's elasticity and provide additional mechanical support (Ouyang et al., 2013; Kirkcaldie and Collins, 2016). The actin cortex generates an overall compressive force around the axon that counteracts axonal tension and will likely affect the development of damage (Fan et al., 2017; García-Grajales et al., 2017). The extracellular matrix and the myelin sheaths around the axon provide additional stability and mechanical support (Goriely et al., 2015b; Weickenmeier et al., 2016). Clearly, further experimental and computational research is needed to qualify and quantify the effects of these structural elements on axonal damage. Second, although the Bell model is widely used for a variety of chemical bonds (Evans and Calderwood, 2007), it is not specific to the tau-microtubule complex. The tau-microtubule interaction is largely unknown and an active field of research (Kadavath et al., 2015; Li et al., 2015; Vemu et al., 2016). An improved understanding of the tau-microtubule binding mechanisms and tau-tau interactions can directly feed into our model and will help improving our model predictions. Third, we assume an S-shaped damage evolution as function of axonal stretch that provides a good representation of the numerical simulation results. However, our current S-shaped curve does not explicitly enforce that a zero damage condition at no loading. In principle, we could use any other damage evolution law to model the stretch- and stretch rate-dependent evolution of damage. Fourth, our method of quantifying axon damage is based on the excessive detachment of crosslinks caused by an externally applied stretch. Within our computational model, however, this axonal damage may recover when the axon is unloaded or when the stretch is held constant and crosslinks reattach to the microtubules. This recovery is non-standard in continuum damage mechanics. Future research should investigate this issue in more detail to improve on the dynamic mechanism assigned to our crosslinks. For example, crosslinks that detach due to an excessive force may not be able or allowed to reattach again, which could be consistent with the experimentally observed increase in tau protein concentration upon axonal damage.

5. Conclusion

The interplay of protein dynamics and physical forces is critical to understand the underlying mechanisms of axonal degradation and brain damage. Here we provides a systematic strategy to relate the discrete dynamic behavior of tau crosslinks on the protein level to the progressive structural degradation on the cellular level to a continuum damage model on the tissue level. Consistent with the definition in nonlinear mechanics, we interpret damage as the gradual stiffness degradation that emerges naturally from a net reduction of crosslinking tau proteins. Motivated by molecular mechanisms, the evolution of damage depends on both the axonal stretch and stretch rate. The only unknown parameter in our model is the characteristic crosslink bond strength, which we vary systematically over two orders of magnitude. Strikingly, for a wide range of stretches, from 1.0 to 1.5, stretch rates, from 0.1 to 10 /s, and bond strengths, from 1 to 100 pN, our model predicts a rather narrow window of critical damage thresholds from 1.01 to 1.30. These values agree well with the experimentally observed axonal damage thresholds reported in the literature. We anticipate that our biophysical model will improve our fundamental understanding of the development and propagation of brain damage across scales and provide useful guidelines to characterize the critical damage level thresholds in response to physical forces.

Author Contributions

RdR performed the simulations. RdR and EK designed the study, analyzed data, and wrote the manuscript.

Funding

This study was supported by the Stanford Graduate Fellowship to RdR and by the Stanford Bio-X IIP seed grant Molecular Mechanisms of Chronic Traumatic Encephalopathy and the NSF grant CMMI 1727268 Understanding Neurodegeneration across the Scales to EK.

Conflict of Interest Statement

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.

References

Ahmadzadeh, H., Smith, D. H., and Shenoy, V. B. (2014). Viscoelasticity of tau proteins leads to strain rate-dependent breaking of microtubules during axonal stretch injury: predictions from a mathematical model. Biophys. J. 106, 1123–1133. doi: 10.1016/j.bpj.2014.01.024

CrossRef Full Text | Google Scholar

Ahmadzadeh, H., Smith, D. H., and Shenoy, V. B. (2015). Mechanical effects of dynamic binding between tau proteins on microtubules during axonal injury. Biophys. J. 109, 2328–2337. doi: 10.1016/j.bpj.2015.09.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Alberts, B., Johnson, A., Lewis, J., Morgan, D., Raff, M., Roberts, K., et al. (2014). Molecular Biology of the Cell, 6th Edn. New York, NY: Garland Science.

Asken, B. M., Sullan, M. J., DeKosky, S. T., Jaffee, M. S., and Bauer, R. M. (2017). Research gaps and controversies in chronic traumatic encephalopathy: a review. JAMA Neurol. 74, 1255–1262. doi: 10.1001/jamaneurol.2017.2396

PubMed Abstract | CrossRef Full Text | Google Scholar

Bain, A. C., and Meaney, D. F. (2000). Tissue-level thresholds for axonal damage in an experimental model of central nervous system white matter injury. J. Biomech. Eng. 122, 615–622. doi: 10.1115/1.1324667

PubMed Abstract | CrossRef Full Text | Google Scholar

Bain, A. C., Raghupathi, R., and Meaney, D. F. (2004). Dynamic stretch correlates to both morphological abnormalities and electrophysiological impairment in a model of traumatic axonal injury. J. Neurotr. 18, 499–511. doi: 10.1089/089771501300227305

PubMed Abstract | CrossRef Full Text | Google Scholar

Beatty, F., and Krishnaswamy, S. (2000). The Mullins effect in equibiaxial deformation. Zeitschrift für angewandte Mathematik Physik ZAMP 51, 984–1015. doi: 10.1007/PL00001533

CrossRef Full Text | Google Scholar

Bell, G. I. (1978). Models for the specific adhesion of cells to cells. Science 200, 618–627. doi: 10.1126/science.347575

PubMed Abstract | CrossRef Full Text | Google Scholar

Bray, D., and Bunge, M. B. (1981). Serial analysis of microtubules in cultured rat sensory axons. J. Neurocytol. 10, 589–605. doi: 10.1007/BF01262592

PubMed Abstract | CrossRef Full Text | Google Scholar

Budday, S., Sommer, G., Hayback, J., Steinmann, P., Holzapfel, G. A., and Kuhl, E. (2017). Rheological characterization of human brain tissue. Acta Biomater. 60, 315–329. doi: 10.1016/j.actbio.2017.06.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Caminiti, R., Carducci, F., Piervincenzi, C., Battaglia-Mayer, A., Confalone, G., Visco-Comandini, F., et al. (2013). Diameter, length, speed, and conduction delay of callosal axons in macaque monkeys and humans: comparing data from histology and magnetic resonance imaging diffusion tractography. J. Neurosci. 33, 14501–14511. doi: 10.1523/JNEUROSCI.0761-13.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, M. C., Chung, P. J., Song, C., Miller, H. P., Kiris, E., Li, Y., et al. (2016). Paclitaxel suppresses tau-mediated microtubule bundling in a concentration-dependent manner. Biochim. Biophys. Acta 1861(1 Pt A), 3456–3463. doi: 10.1016/j.bbagen.2016.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Chung, P. J., Choi, M. C., Miller, H. P., Feinstein, H. E., Raviv, U., Li, Y., et al. (2015). Direct force measurements reveal that protein tau confers short-range attractions and isoform-dependent steric stabilization to microtubules. Proc. Natl. Acad. Sci. USA 112, E6416–E6425. doi: 10.1073/pnas.1513172112

PubMed Abstract | CrossRef Full Text | Google Scholar

Cloots, R. J. H., van Dommelen, J. A. W., Kleiven, S., and Geers, M. G. D. (2013). Multi-scale mechanics of traumatic brain injury: predicting axonal strains from head loads. Biomech. Model. Mechanobiol. 12, 137–150. doi: 10.1007/s10237-012-0387-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Coles, C. H., and Bradke, F. (2015). Coordinating neuronal actin–microtubule dynamics. Curr. Biol. 25, R677–R691. doi: 10.1016/j.cub.2015.06.020

PubMed Abstract | CrossRef Full Text | Google Scholar

de Borst, R., Sluys, L. J., Muhlhaus, H. B., and Pamin, J. (1993). Fundamental issues in finite element analyses of localization of deformation. Eng. Comput. Int. J. Comput. Aided Eng. 10, 99–121. doi: 10.1108/eb023897

CrossRef Full Text | Google Scholar

de Rooij, R., and Kuhl, E. (2016). Constitutive modeling of brain tissue: current perspectives. Appl. Mech. Rev. 68, 1–16. doi: 10.1115/1.4032436

CrossRef Full Text | Google Scholar

de Rooij, R., and Kuhl, E. (2018). Microtubule polymerization and cross-link dynamics explain axonal stiffness and damage. Biophys. J. 114, 201–212. doi: 10.1016/j.bpj.2017.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

de Rooij, R., Miller, K. E., and Kuhl, E. (2016). Modeling molecular mechanisms in the axon. Comput. Mech. 59, 523–537. doi: 10.1007/s00466-016-1359-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Debanne, D., Campanac, E., Bialowas, A., Carlier, E., and Alcaraz, G. (2011). Axon physiology. Physiol. Rev. 91, 555–602. doi: 10.1152/physrev.00048.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Dogterom, M., and Yurke, B. (1997). Measurement of the force-velocity relation for growing microtubules. Science 278, 856–860. doi: 10.1126/science.278.5339.856

PubMed Abstract | CrossRef Full Text | Google Scholar

Eisenberg, D. S., and Sawaya, M. R. (2017). Neurodegeneration: taming tangled tau. Nature 547, 170–171. doi: 10.1038/nature23094

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans, E., and Ritchie, K. (1997). Dynamic strength of molecular adhesion bonds. Biophys. J. 72, 1541–1555. doi: 10.1016/S0006-3495(97)78802-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans, E. A., and Calderwood, D. A. (2007). Forces and bond dynamics in cell adhesion. Science 316, 1148–1153. doi: 10.1126/science.1137592

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, A., Tofangchi, A., Kandel, M., Popescu, G., and Saif, T. (2017). Coupled circumferential and axial tension driven by actin and myosin influences in vivo axon diameter. Sci. Rep. 7:14188. doi: 10.1038/s41598-017-13830-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Fitzpatrick, A. W. P., Falcon, B., He, S., Murzin, A. G., Murshudov, G., Garringer, H. J., et al. (2017). Cryo-EM structures of tau filaments from Alzheimer's disease. Nature 56:343. doi: 10.1038/nature23002

CrossRef Full Text | Google Scholar

García-Grajales, J. A., Jérusalem, A., and Goriely, A. (2017). Continuum mechanical modeling of axonal growth. Comput. Methods Appl. Mech. Eng. 314, 147–163. doi: 10.1016/j.cma.2016.07.032

CrossRef Full Text

Geers, M. G. D., Brekelmans, W. A. M., and de Borst, R. (1994). “Viscous regularization of strain-localisation for damaging materials,” in DIANA Computational Mechanics ‘94 (Dordrecht: Springer), 127–138.

Google Scholar

Giordano, C., Zappalà, S., and Kleiven, S. (2017). Anisotropic finite element models for brain injury prediction: the sensitivity of axonal strain to white matter tract inter-subject variability. Biomech. Model. Mechanobiol. 16, 1269–1293. doi: 10.1007/s10237-017-0887-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Gittes, F., Mickey, B., Nettleton, J., and Howard, J. (1993). Flexural rigidity of microtubules and actin filaments measured from thermal fluctuations in shape. J. Cell Biol. 120, 923–934. doi: 10.1083/jcb.120.4.923

PubMed Abstract | CrossRef Full Text | Google Scholar

Goriely, A., Budday, S., and Kuhl, E. (2015a). Neuromechanics: from neurons to brain. Adv. Appl. Mech. 48, 79–139. doi: 10.1016/bs.aams.2015.10.002

CrossRef Full Text | Google Scholar

Goriely, A., Geers, M. G. D., Holzapfel, G. A., Jayamohan, J., Jerusalem, A., Sivaloganathan, S., et al. (2015b). Mechanics of the brain: perspectives, challenges, and opportunities. Biomech. Model. Mechanobiol. 14, 931–965. doi: 10.1007/s10237-015-0662-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Greenwald, R. M., Gwin, J. T., Chu, J. J., and Crisco, J. J. (2008). Head impact severity measures for evaluating mild traumatic brain injury risk exposure. Neurosurgery 62, 789–798. doi: 10.1227/01.neu.0000318162.67472.ad

PubMed Abstract | CrossRef Full Text | Google Scholar

Haak, R. A., Kleinhans, F. W., and Ochs, S. (1976). The viscosity of mammalian nerve axoplasm measured by electron spin resonance. J. Physiol. 263, 115–137. doi: 10.1113/jphysiol.1976.sp011624

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirokawa, N. (1982). Cross-linker system between neurofilaments, microtubules and membranous organelles in frog axons revealed by the quick-freeze, deep-etching method. J. Cell Biol. 94, 129–142. doi: 10.1083/jcb.94.1.129

PubMed Abstract | CrossRef Full Text | Google Scholar

Howard, J. (2001). Mechanics of Motor Proteins and the Cytoskeleton. Sunderland: Sinauer Associates, Inc.

Google Scholar

Hyder, A. A., Wunderlich, C. A., Puvanachandra, P., Gururaj, G., and Kobusingye, O. C. (2007). The impact of traumatic brain injuries: a global perspective. NeuroRehabilitation 22, 341–353.

PubMed Abstract | Google Scholar

Hyland, C., Mertz, A. F., Forscher, P., and Dufresne, E. (2014). Dynamic peripheral traction forces balance stable neurite tension in regenerating Aplysia bag cell neurons. Sci. Rep. 4:4961. doi: 10.1038/srep04961

PubMed Abstract | CrossRef Full Text | Google Scholar

Igaev, M., Janning, D., Sündermann, F., Niewidok, B., Brandt, R., and Junge, W. (2014). A refined reaction-diffusion model of tau-microtubule dynamics and its application in FDAP analysis. Biophys. J. 107, 2567–2578. doi: 10.1016/j.bpj.2014.09.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Jakobs, M., Franze, K., and Zemel, A. (2015). Force generation by molecular-motor-powered microtubule bundles; implications for neuronal polarization and growth. Front. Cell. Neurosci. 9:441. doi: 10.3389/fncel.2015.00441

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, S., Zhao, W., Li, Z., and McAllister, T. W. (2014). Head impact accelerations for brain strain-related responses in contact sports: a model-based investigation. Biomech. Model. Mechanobiol. 13, 1121–1136. doi: 10.1007/s10237-014-0562-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, V. E., Stewart, W., and Smith, D. H. (2013). Axonal pathology in traumatic brain injury. Exp. Neurol. 246, 35–43. doi: 10.1016/j.expneurol.2012.01.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Kachanov, L. M. (1986). Introduction to Continuum Damage Mechanics. Dordrecht: Martinus Nijhoff Publishers.

Google Scholar

Kadavath, H., Jaremko, M., Jaremko, Ł., Biernat, J., Mandelkow, E., and Zweckstetter, M. (2015). Folding of the tau protein on microtubules. Angew. Chem. Int. Edn. 54, 10347–10351. doi: 10.1002/anie.201501714

PubMed Abstract | CrossRef Full Text | Google Scholar

Kirkcaldie, M. T. K., and Collins, J. M. (2016). The axon as a physical structure in health and acute trauma. J. Chem. Neuroanat. 76, 9–18. doi: 10.1016/j.jchemneu.2016.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Kleiven, S., and von Holst, H. (2002). Consequences of head size following trauma to the human head. J. Biomech. 35, 153–160. doi: 10.1016/S0021-9290(01)00202-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Krieg, M., Stühmer, J., Cueva, J. G., Fetter, R., Spilker, K., Cremers, D., et al. (2017). Genetic defects in β-spectrin and tau sensitize C. elegans axons to movement-induced damage via torque-tension coupling. eLife 6:e20172. doi: 10.7554/eLife.20172

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuhl, E., Ramm, E., and de Borst, R. (2000). An anisotropic gradient damage model for quasi-brittle materials. Comput. Methods Appl. Mech. Eng. 183, 87–103. doi: 10.1016/S0045-7825(99)00213-3

CrossRef Full Text | Google Scholar

Kuo, C., Wu, L. C., Ye, P. P., Laksari, K., Camarillo, D. B., and Kuhl, E. (2017). Pilot findings of brain displacements and deformations during roller coaster rides. J. Neurotr. 34, 3198–3205. doi: 10.1089/neu.2016.4893

PubMed Abstract | CrossRef Full Text | Google Scholar

Lazarus, C., Soheilypour, M., and Mofrad, M. R. K. (2015). Torsional behavior of axonal microtubule bundles. Biophys. J. 109, 231–239. doi: 10.1016/j.bpj.2015.06.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Lemaitre, J. (1992). A Course on Damage Mechanics. Berlin; Heidelberg; New York: Springer.

Google Scholar

Li, W. (2016). Damage models for soft tissues: a survey. J. Med. Biol. Eng. 36, 285–307. doi: 10.1007/s40846-016-0132-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X.-H., Culver, J. A., and Rhoades, E. (2015). Tau binds to multiple tubulin dimers with helical structure. J. Am. Chem. Soc. 137, 9218–9221. doi: 10.1021/jacs.5b04561

PubMed Abstract | CrossRef Full Text | Google Scholar

Mallik, R., Carter, B. C., Lex, S. A., King, S. J., and Gross, S. P. (2004). Cytoplasmic dynein functions as a gear in response to load. Nature 427, 649–652. doi: 10.1038/nature02293

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, H., Zhang, L., Jiang, B., Genthikatti, V. V., Jin, X., Zhu, F., et al. (2013). Development of a finite element human head model partially validated with thirty five experimental cases. J. Biomech. Eng. 135:111002. doi: 10.1115/1.4025101

PubMed Abstract | CrossRef Full Text | Google Scholar

Margulies, S. S., and Thibault, L. E. (1992). A proposed tolerance criterion for diffuse axonal injury in man. J. Biomech. 25, 917–923. doi: 10.1016/0021-9290(92)90231-O

PubMed Abstract | CrossRef Full Text | Google Scholar

Marini, G., Maier, A., Reeps, C., Eckstein, H. H., Wall, W. A., and Gee, M. W. (2012). A continuum description of the damage process in the arterial wall of abdominal aortic aneurysms. Int. J. Numer. Methods Biomed. Eng. 28, 87–99. doi: 10.1002/cnm.1472

PubMed Abstract | CrossRef Full Text | Google Scholar

Mez, J., Daneshvar, D. H., Kiernan, P. T., Abdolmohammadi, B., Alvarez, V. E., Huber, B. R., et al. (2017). Clinicopathological evaluation of Chronic Traumatic Encephalopathy in players of American football. JAMA 318, 360–370. doi: 10.1001/jama.2017.8334

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Toole, M., Lamoureux, P., and Miller, K. E. (2008). A physical model of axonal elongation: force, viscosity, and adhesions govern the mode of outgrowth. Biophys. J. 94, 2610–2620. doi: 10.1529/biophysj.107.117424

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Toole, M., Lamoureux, P., and Miller, K. E. (2015). Measurement of subcellular force generation in neurons. Biophys. J. 108, 1027–1037. doi: 10.1016/j.bpj.2015.01.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Ouyang, H., Nauman, E., and Shi, R. (2013). Contribution of cytoskeletal elements to the axonal mechanical properties. J. Biol. Eng. 7, 21.1–21.8. doi: 10.1186/1754-1611-7-21

PubMed Abstract | CrossRef Full Text | Google Scholar

Peerlings, R. H. J., Geers, M. G. D., de Borst, R., and Brekelmans, W. A. M. (2001). A critical comparison of nonlocal and gradient-enhanced softening continua. Int. J. Solids Struct. 38, 7723–7746. doi: 10.1016/S0020-7683(01)00087-7

CrossRef Full Text | Google Scholar

Pereira, L. F., Weerheijm, J., and Sluys, L. J. (2017). A new effective rate dependent damage model for dynamic tensile failure of concrete. Eng. Fract. Mech. 176, 281–299. doi: 10.1016/j.engfracmech.2017.03.048

CrossRef Full Text | Google Scholar

Peter, S. J., and Mofrad, M. R. K. (2012). Computational modeling of axonal microtubule bundles under tension. Biophys. J. 102, 749–757. doi: 10.1016/j.bpj.2011.11.4024

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajagopalan, J., Tofangchi, A., and A Saif, M. T. (2010). Drosophila neurons actively regulate axonal tension in vivo. Biophys. J. 99, 3208–3215. doi: 10.1016/j.bpj.2010.09.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Recho, P., Jérusalem, A., and Goriely, A. (2016). Growth, collapse, and stalling in a mechanical model for neurite motility. Phys. Rev. E 93:032410. doi: 10.1103/PhysRevE.93.032410

CrossRef Full Text | Google Scholar

Singh, A., Lu, Y., Chen, C., Kallakuri, S., and Cavanaugh, J. M. (2016). A new model of traumatic axonal injury to determine the effects of strain and displacement rates. Stapp Car Crash J. 50, 601–623.

PubMed Abstract | Google Scholar

Smith, D. H., and Meaney, D. F. (2016). Axonal damage in traumatic brain injury. Neuroscientist 6, 483–495. doi: 10.1177/107385840000600611

CrossRef Full Text | Google Scholar

Suresh, S. (2007). Biomechanics and biophysics of cancer cells. Acta Biomater. 3, 413–438. doi: 10.1016/j.actbio.2007.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Suter, D. M., and Miller, K. E. (2011). The emerging role of forces in axonal elongation. Prog. Neurobiol. 94, 91–101. doi: 10.1016/j.pneurobio.2011.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Takhounts, E. G., Eppinger, R. H., Campbell, J. Q., Tannous, R. E., Power, E. D., and Shook, L. S. (2003). On the development of the SIMon finite element head model. Stapp Car Crash J. 47, 107–133.

PubMed Abstract | Google Scholar

Tang-Schomer, M. D., Johnson, V. E., Baas, P. W., Stewart, W., and Smith, D. H. (2012). Partial interruption of axonal transport due to microtubule breakage accounts for the formation of periodic varicosities after traumatic axonal injury. Exp. Neurol. 233, 364–372. doi: 10.1016/j.expneurol.2011.10.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang-Schomer, M. D., Patel, A. R., Baas, P. W., and Smith, D. H. (2010). Mechanical breaking of microtubules in axons during dynamic stretch injury underlies delayed elasticity, microtubule disassembly, and axon degeneration. FASEB J. 24, 1401–1410. doi: 10.1096/fj.09-142844

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, C. A. (2017). Traumatic brain injury–related emergency department visits, hospitalizations, and deaths — United States, 2007 and 2013. MMWR Surveil. Summ. 66, 1–16. doi: 10.15585/mmwr.ss6609a1

PubMed Abstract | CrossRef Full Text | Google Scholar

Tofangchi, A., Fan, A., and Saif, M. T. A. (2016). Mechanism of axonal contractility in embryonic drosophila motor neurons in vivo. Biophys. J. 111, 1519–1527. doi: 10.1016/j.bpj.2016.08.024

PubMed Abstract | CrossRef Full Text | Google Scholar

van den Bedem, H., and Kuhl, E. (2015). Tau-ism: the yin and yang of microtubule sliding, detachment, and rupture. Biophys. J. 109, 2215–2217. doi: 10.1016/j.bpj.2015.10.020

PubMed Abstract | CrossRef Full Text | Google Scholar

van den Bedem, H., and Kuhl, E. (2017). Molecular mechanisms of chronic traumatic encephalopathy. Curr. Opin. Biomed. Eng. 1, 23–30. doi: 10.1016/j.cobme.2017.02.003

CrossRef Full Text | Google Scholar

Vemu, A., Atherton, J., Spector, J. O., Szyk, A., Moores, C. A., and Roll-Mecak, A. (2016). Structure and dynamics of single-isoform recombinant neuronal human tubulin. J. Biol. Chem. 291, 12907–12915. doi: 10.1074/jbc.C116.731133

PubMed Abstract | CrossRef Full Text | Google Scholar

Verhulst, P. F. (1845). Recherches mathématiques sur la loi d'accroissement de la population. Nouveaux Mémoires de l'Académie Royale Sciences Belles-Lettres Bruxelles 18, 1–42.

Google Scholar

Wegmann, S., Scholer, J., Bippes, C. A., Mandelkow, E., and Muller, D. J. (2011). Competing interactions stabilize pro- and anti-aggregant conformations of human tau. J. Biol. Chem. 286, 20512–20524. doi: 10.1074/jbc.M111.237875

PubMed Abstract | CrossRef Full Text | Google Scholar

Weickenmeier, J., Butler, C. A. M., Young, P. G., Goriely, A., and Kuhl, E. (2017). The mechanics of decompressive craniectomy: personalized simulations. Comput. Methods Appl. Mech. Eng. 314, 180–195. doi: 10.1016/j.cma.2016.08.011

CrossRef Full Text | Google Scholar

Weickenmeier, J., de Rooij, R., Budday, S., Steinmann, P., Ovaert, T. C., and Kuhl, E. (2016). Brain stiffness increases with myelin content. Acta Biomater. 42, 265–272. doi: 10.1016/j.actbio.2016.07.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Weisbecker, H., Pierce, D. M., Regitnig, P., and Holzapfel, G. A. (2012). Layer-specific damage experiments and modeling of human thoracic and abdominal aortas with non-atherosclerotic intimal thickening. J. Mech. Behav. Biomed. Mater. 12, 93–106. doi: 10.1016/j.jmbbm.2012.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Woerman, A. L., Aoyagi, A., Patel, S., Kazmi, S. A., Lobach, I., Grinberg, L. T., et al. (2016). Tau prions from Alzheimer's disease and chronic traumatic encephalopathy patients propagate in cultured cells. Proc. Nat. Acad. Sci. 113, E8187–E8196. doi: 10.1073/pnas.1616344113

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, W., and Baas, P. W. (1994). Changes in microtubule number and length during axon differentiation. J. Neurosci. 14(5 Pt 1), 2818–2829. doi: 10.1523/JNEUROSCI.14-05-02818.1994

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuen, T. J., Browne, K. D., Iwata, A., and Smith, D. H. (2009). Sodium channelopathy induced by mild axonal trauma worsens outcome after a repeat injury. J. Neurosci. Res. 87, 3620–3625. doi: 10.1002/jnr.22161

PubMed Abstract | CrossRef Full Text | Google Scholar