Mathematical modelling of microtubule-tau protein transients: insights into the superior mechanical behavior of axon

Axons with staggered microtubules cross-linked by tau protein possess a remarkable mechanical balance of high specific stiffness and toughness. Owing to their viscoelastic nature, axons exhibit stress rate-dependent mechanical behavior, which is relevant to their selective vulnerability to damage in traumatic brain injury. A Kelvin-Voigt viscoelastic shear lag model is developed to elucidate the mechanical responses of axons under transient tensile force. Analytical closed-form expressions are derived to characterize the relative sliding, stress transfer and failure mechanism between microtubule and tau protein while the axon is stretched transiently. The results from the theoretical solutions elucidate how the MT-tau interface length and stress rate affect the mechanical responses of axon. It is found that axonal failure mechanism may be different under different loading conditions. Long microtubules are more vulnerable to rupture at high stress rate, yet short microtubules are likely to detach from microtubule bundles under large deformations. In the view of multi-level failure of axon, it is illustrated how the vulnerable axons protect themselves from overall damage, and how the axon can simultaneously achieve an outstanding mechanical balance of high specific stiffness and toughness.


Introduction
Staggered biocomposites are found to exhibit a remarkable mechanical balance of fracture toughness versus stiffness and strength [1,2]. In consideration of the structure of staggered biocomposite and the vulnerability of axon, it may possess the innate ability of mechanical balance of high stiffness and high toughness, and self-protecting over the course of evolution; hence the axon is taken as an example for investigation in this study. Traumatic brain injury (TBI) is primarily resulted from an impact to the brain tissue from an external mechanical force including contact force, penetration of a projectile, or inertia forces induced by rapid 2 acceleration and deceleration [3]. In recent years, TBI has become a major public health issue, killing 15 per 100,000 individuals in European countries [4] and affecting over 2.5 million people in the USA [5] per year.
Plenty of in vitro researches indicated that the mechanical rupture of microtubules (MTs) results in traumatic brain injury, commonly referred to as "concussion" [6,7]. As depicted in Fig. 1a, a single long slender axon and multiple shorter dendrites are responsible for transmitting and receiving signals, respectively [8]. Histologically, axonal cytoskeleton consists of MT bundles connected by microtubule-associated protein (MAP) tau, as illustrated in Fig. 1b [8,9]. Parallel array of MTs, acting as reinforcing components with a main function of transporting electrochemical cargo, are the stiffest structural element within the axon [10][11][12][13].
Owing to their viscoelastic nature, axons exhibit strain or stress rate-dependent rupture [7]. In order to characterize the micromechanical breaking of axonal ultrastructure, axons were uniaxially stretched transiently by using transmission electron microscopy (TEM) under different strain rates [6,7].
Experimental techniques are tricky to capture the immediate mechanical behavior of axonal MT bundles in tension. It is reported that the axon behaved like an elastic element while it is subjected to sudden changes in force [14]. Indeed, viscoelasticity may be more appropriate to capture the major biophysical behavior of axon and associated MTs in such conditions.
In recent years, there is an increasing interest in modelling the mechanical response of axon. However, little is known about which component of the axonal cytoskeleton might break under transient mechanical loading and large deformation. Tolomeo et al. [15] reported that cross-linking proteins provided shear resistance between MTs. Peter et al. [16] used a discrete bead-spring model to simulate the biomechanical behavior of axonal microtubule bundle under uniaxial tension, and they assumed that both MT and tau protein are linear elastic materials. Shamloo et al. [17] considered MTs as a large number of discrete masses and modelled tau protein as Kelvin-Voigt element so that viscoelastic model can be employed to simulate the transient response of axonal MTs under sudden forces. de Rooij et al.
modelled the time-dependent behaviour of brain tissue by adopting an analytical hyperelastic approach [18], and developed a computational mechanics model to elucidate cellular-level characteristics of MT bundles [19]. Combined with experimental results and finite element analysis (FEA), viscosity has been identified as the most critical factor causing the mechanical vulnerability of axon under transient loading [10].
Within the axonal cytoskeleton, tau protein plays a vital role in the assembly of individual MTs [20], and the breaking of tau protein may be relevant to some of neurodegenerative diseases, such as Alzheimer's disease [21,22] and Parkinsonism [23]. Different viscoelastic models consisting of linear elastic MT and viscoelastic tau protein, a combination of elastic springs and viscous dashpots, were postulated to elucidate the macroscopic viscoelastic behavior of the axon by some researchers [14,24,25]. Ahmadzadeh et al. [10] verified that Kelvin viscoelastic model can efficiently characterize the viscoelastic mechanical behavior of axon, but only MT failure is investigated in their model without consideration of tau protein breaking. However, once tau protein axial strain exceeds its ultimate strain, MTs will detach from the MT bundle [26,27]. What's more, tau protein breaking is found to be associated with some of neurodegenerative diseases. Therefore, it's significant to develop a 3 computational model comprehensively characterizing the axonal failure evidence of MT rupture and tau protein breaking.
Though many researches revealed the mechanical behavior of tau protein and MTs [10,14,24,25], the molecular mechanism of axonal failure remains poorly understood. In Ref. [28]  3. How does the axon simultaneously achieve an outstanding mechanical balance of high specific stiffness and toughness?

Fundamental formulas and theoretical derivation a b c
In 1952, SLM was firstly proposed by Cox [30], which has now been modified and developed for biomechanical applications [31][32][33][34][35][36][37][38]. SLM can be employed to characterize stress transfer between MTs and tau protein under transient tensile conditions. As shown in Fig. 2a, the cytoskeleton of MT is hollow cylindrical structures with outer radius RO and inner radius RI. As the elastic modulus of MTs is much larger than that of tau protein, MTs are assumed to follow linear elastic constitutive relation, and its Young's modulus is denoted by EM.
Before theoretical derivations, the assumptions are adopted in the analysis presented in this study as follows: (i) The stiffest discontinuous MTs are connected by tau protein, leading to deformation and deformation rate primarily occur in tau protein, thus, all the MTs are assumed to be linear elastic materials, and tau protein follows viscoelastic constitutive law.
(ii) The MT-tau interface length L is over two orders of magnitude larger than the distance between adjacent MTs so that the deformation in axon is essentially along x-direction.
(iii) All the MTs are parallel array along x-direction, and the MT-tau interface length L is nearly half of their entire length (see Fig. 2a). Recently, researchers found that axons have viscoelastic material behavior [5,17,39,40]. Based on the experiment [41], a Kelvin-Voigt viscoelastic model containing a spring with stiffness K in parallel with a dashpot with viscosity µ (Fig. 2d) is employed for tau protein in order to predict the biomechanical behavior of axon subjected to transient tensile force where η=µ/K is tau protein dashpot timescale. A Cartesian coordinate system is established at the midpoint of an entire length of MT so that the longitudinal direction of MTs is along x-axis (see Fig.   5 2c). In the existing study of axon [10], the axon is simplified as a staggered arrangement of discontinuous MTs connected by tau proteins, which is represented by a periodic unit cell that is half the entire length of MT. When the unit cell is stretched by a tensile force, stress will be transferred from one MT to the adjacent MTs through tau protein, which leads to relative sliding between MTs. Tau protein elongation δ and the relative sliding between adjacent MTs [u1(x,t)-u2(x,t)] follow the geometrical relationship as follows where β is the angle between the tau protein and the MT; u1 and u2 are displacements of upper and lower MTs, respectively.
Tensile forces F(x,t) resulting from the elongation of tau protein can be converted to shear stress τ(x,t) acting over the entire MT circumference, which can be written as where dT is the center-to-center distance between neighboring tau protein. The number of adjacent MTs is α=6 based on the cross-sectional electron micrograph of MT bundle ( Fig. 1c and Fig. 2b).
According to the unit cell (see Fig. 2c) showing the two adjacent MTs, shear stress τ(x,t) acts over the MT circumference leading to the following equilibrium equations where RO and RI are the outer radius and inner radius of the MTs, respectively. σ1 and σ2 are the normal stresses of the linear elastic upper and lower MTs, which can be expressed as where EM is the Young's modulus of MTs. ε1 and ε2 are the axial strain of the upper and lower MTs, respectively.
While MT has a negligible mass, at any point x along MTs longitude direction, the tensile force P(t) should maintain mechanical balance with the normal stresses acting on the MT cross-section, thus we can obtain where and combining Eqs. (1)-(8), we can get the governing equation for tau protein elongation 2 2 The upper MT normal stress σ1 along the x-direction can be written as Based on Eq. (10), boundary and initial conditions are alternatively expressed as where LC is the characteristic length over which the stress is transferred between MT and tau protein.
Tau protein elongation can be derived from solving the nonhomogeneous partial differential equation (Eq. (9)) with nonhomogeneous boundary conditions (Eq. (11)) and initial condition (Eq. (12)) It can be seen from the above equations that, the values of C0 and Cn depend on the initial condition. If there are no elongation and no force in the beginning, we can get C0=Cn=0. D0 and Dn are determined Obviously, MT axial strain can be derived from Eq. (6). Tau protein axial strain can be derived from the definition of axial strain εT=δsinβ/dM, where dM is the surface-to-surface distance between adjacent 7 MTs.

Failure criteria and universal applicability of the solutions
In general, the mechanical properties of strength and toughness are mutually exclusive [42]; yet these properties are a significant requirement for most structural materials. Strength is invariably a stress representing a material's resistance to non-recoverable deformation [43]. When the transient tensile force continuously exerts at the end of the MT, the main component of either MT or tau protein will reach the ultimate strain. It was reported that MTs are shown to be ruptured under 50% of axial strain by several researchers [44][45][46]. With increasing strain, tau-tau bonds begin to break in exceed of approximately 40% [41].

Related experiments and parameters
It can be noted from Eq. (9) that, the micromechanical interaction behavior of axon is determined by (LC/ζ) 2  Experimental data from numerous researches are applied to assign values to the parameters in viscoelastic SLM [9,16,17,41,[44][45][46][47][48][49]. The length of MT is about 1-10 μm, and its outer and inner diameters are 25 nm and 14 nm, respectively [9]. Axonal MTs are hexagonally aligned with surface-to-surface spacing of approximately 20 nm [16,17,29]. In general, the Young's modulus of MT could be captured by the slope of its stress-strain curve. The ultimate tensile strain of MT can be defined as the axial strain when the MT ruptures. The Young's modulus of MTs was found to be 1.51 ± 0.19 GPa [9], which agreed well with experimental data reported by Suresh [50]. The rupture tensile strain can be set at 50% according to experimental measurements of the MT rupture strain performed by Janmey et al. [44], consistent with Suresh's experimental observation [50].
The geometrical and material parameters of tau protein are also of great significance. About 95% of tau proteins are formed by paired helical filaments with a radius of 4-10 nm [47]. Based on statistical analysis, Hirokawa et al. [48] reported that tau protein center-to-center spacing fell within the typical range, 20-40 nm. The Young's modulus of tau protein is typically about 5 MPa [16,49], which gives the tau cross-link axial spring constant of 33.3 pN/nm as calculated from K=ETπR 2 T sinβ/dM. According to the experimental data [41], Ahmadzadeh et al. [10] estimated the viscoelastic parameter η=0.35 s by using Bell's equation [51]. Additionally, interactions in hTau40 broke at maximum force 300 pN observed by Wegmann et al. [41], and therefore tau protein breaking tensile strain is estimated in exceed of 40%.
In the current work, the geometrical and material parameters are gained originally from experiment data or estimations, as listed in Table 1. The angle β between the tau protein and the MT is assumed to be 60°. According to these parameters, we can obtain the characteristic length from Eq. (8), LC=0.390 μm.    continuously stretched at a high stress rate, tau protein axial strain is extremely large at both ends (x=0 or L), but is nearly zero in the region 0.1<x/L<0.9 ( Fig. 5a and 5c). The axial strain distribution of the upper MT and lower MT surrounding the midpoint (0.1<x/L<0.9) remain constant, which equals to half of the maximum strain ( Fig. 5b and 5d). It can be noted from Fig. 5b and 5d that, evidence of MT a b

Sliding-rupture regime
rupture is more sensitive to stress rate rather than tensile stress, consistent with the experimental observation in the region of axonal swelling by Tang-Schomer et al. [7]. MTs are considered as the tracks for chemical cargo transport, and they will rupture under high stress rate (  ≥92.5EM Pa/s) found in TBI, leading to the transport interruption [6,7]. This mathematical modelling can help us to understand the mechanism of TBI, and predict the symptom appearance of TBI when the brain is subjected to sudden external mechanical forces.  Fig. 2c, respectively).

Sliding-detachment regime
The second type of axonal failure is caused by tau protein breaking, leading to MT detachment from MT bundle. Tau protein and MT axial strain are shown in Fig. 6. When applied stress rate are much smaller than 92.5EM Pa/s, e.g.  =10EM Pa/s or 0.1EM Pa/s, axonal failure follows S-D regime. It can be seen from Fig. 6 that, as the stretched loading continuously increases, the breaking of bonds between the paired tau proteins starts from both ends and propagates sequentially toward the midpoint of the MT-tau interface, as observed in the experiment [16]. According to the MT detachment time, MT-tau interface breaks faster and faster, as the applied loading increases continuously.

Sliding-detachment-rupture regime
A complicated failure process of axon, containing tau protein breaking and MT rupture, is illustrated in Fig. 7. When applied stress rate is slightly smaller than 92.5EM Pa/s, e.g.  =50EM Pa/s or 75EM Pa/s, shear stress will be concentrated at both ends of the MT-tau interface even though the relative sliding is small. Obviously, tau protein closed to both ends breaks first, and then MTs rupture not until all the connected tau proteins come to failure, as shown in Fig. 7.
Both MT rupture and tau protein breaking will occur in different MT-tau interface lengths under transient loading conditions [3,16]. However, to the best of our knowledge, there is no experimentally

Effect of the stress rate on axonal failure process
In this section, we investigate the strain distribution at different stress rates by comparing Figs. 5-7.
While axon is stretched under  >50EM Pa/s, axial strain of tau protein at both ends (x=0 and L) is extremely large and tau protein axial strain in the region 0.1<x/L<0.9 is close to zero; the MT axial strain distribution surrounding the midpoint (x=L/2) remains constant. While the applied stress rate is very small, the tau protein axial strain distribution surrounding the midpoint (x=L/2) is not close to zero, and MT axial strain surrounding the midpoint is variable, as shown in Fig. 6a. Thus, we can conclude that during the loading process MT axial strain along the interface length is nearly unchanged, and tau-tau bonds break simultaneously in the end if the applied stress rate is small enough (quasi-static).

a b c d
In order to clarify the relationship between the applied stress rate and mechanical response of axons containing Kelvin-Voigt viscoelastic tau protein, the time and spatial coordinate system should be rescaled as Tt   , X x L  , in such a way that variables in Eq. (9) are normalized as According to the rescale variables, Eq. (9) can be converted to Eq. (16) implies that tau protein with high viscosity has similar relative sliding with high applied stress rate when all other parameters remain the same. On the contrary, tau protein with high spring constant has an opposite tendency. Interestingly, axon bundle can withstand strains of over 100% and recover to its original configuration without evidence of damage [7]. However, axons appeared to failure when subjected to a sudden change in tensile force [52]. Yuen et al. [53] reported that MT rupture begins at axon strain 5%, when axon is stretched under strain rate 50 s -1 . This viscoelastic shear lag model can explain why the axonal injury is sensitive to the stress rate.

Effect of the MT-tau interface length on stress transfer
Researchers found that the interface length affects the stress transferring and fracture propagation between fiber and matrix in the composites and biocomposites [54,55]. Recently, finite element models were established to explore how MT-tau interface length affects the mechanical behavior of axon [10]. In this study, SLM including viscoelastic tau protein is constructed in order to illustrate the effects of MT-tau interface length on relative sliding and stress transfer under transient tensile force, as shown in Fig. 8.
While MTs are stretched by a force with a low stress rate  =0.001EM Pa/s, tau protein axial strain along short MT-tau interface length (L=LC) is very large (Fig. 8c). Whereas if the applied stress rate is 10EM Pa/s, tau protein axial strains in both ends are identical, which indicates that the initiation of tau protein breaking starts simultaneously, and overall tau proteins along the short interface length come to failure in a shorter period of time (Fig. 8a). Normal stress in the MT will be concentrated at the both ends when the stress rate is high and MT-tau interface is long enough, as illustrated in Fig. 8b and 8d.
Histologically, long and short MTs are parallel aligned in the axon bundle. Bass et al. [56] found that different lengths of MT have different mechanical functions. Long MT is responsible for carrying loads, which results in a higher resistance to the sliding [57], and implies that long MTs are more susceptible to rupture. Yet, tau proteins in the short MT-tau interface are responsible for sliding, which may explain the significant elongation of axon following traumatic stretch. Our model prediction is consistent with  force-displacement relation of axon, obtained from the experiments, can be found in Refs. [14,58].

Effective stiffness of axon
A series of in vivo tensile tests on axon were recently conducted to measure the mechanical properties c d of axon [14,[59][60][61]. Impressively, Rajagopalan et al. [14] reported that the axons behave like elastic springs under suddenly applied force condition. The effective stiffness of the axons, represented by the slope of the force-deformation curve, takes values in the range 0.2-1.2 nN/μm. The force-deformation curve of axonal bundle under tensile stress is plotted in Fig. 9. It shows that the slope of force-elongation curve of our current model agrees well with that of experimental data, which implies that our prediction for the effective stiffness of the axonal bundles is in excellent agreement with the existing experiment.

Discussion
Diffuse axonal injury (DAI) is one of the most common types of TBI [53]. Unfortunately, today's imaging methods are unable to detect the microscopic diffuse injuries to the axon [62][63][64]. Moreover, tension of axon bundle is such a tricky experiment, which has been prompting scientists to resort to computational modelling. Computational modelling is considered as an efficient tool to predict the mechanical response of axons for substantially improving contemporary understanding of the underlying pathology and molecular mechanisms of axonal injury.  Fig. 10. Therefore, the vulnerable axon can achieve high stiffness and toughness, protecting itself from overall damage. Fig. 10. Different types of deformation and damage evolution in the axon after TBI due to various MT-tau interface lengths and stress rates. Tau protein serves as a molecular switch between MT sliding, rupture, detachment and detachment-rupture.
In our previous study [65], we have explained that the fundamental load-carrying elements of regularly staggered biocomposite structures are discontinuous fibers, instead of continuous long fibers along the entire tissue. The discontinuous fibers are more capable to protect themselves from overall damage. As for axons, MTs and tau proteins will come to failure at different locations in the axon bundle, rather than break in a single section, and therefore an overall mechanical stability and high toughness could be achieved. This may be a key reason why randomly staggered biocomposite structures are observed in the axon.
The mechanical properties of engineering fiber reinforced composites depend to a great extent on the bond length. In contrast to the axon, high-performance fiber reinforced composites typically use continuous fibers, thus achieving high stiffness and strength but presenting limited toughness and ductility. From the above analysis, we know that the use of discontinuous fibers could potentially improve the ductility and fracture toughness [66,67]. Longer interface length is not conducive to stress transfer, but discontinuous fibers could help dissipate energy and protect biocomposites from overall damage. Additionally, crack bridging by discontinuous fibers can make brittle materials tougher by transferring stress from the crack tip to the ductile matrix.
Despite the ultrastructure of the axon has the advantage of dissipating energy, MTs are vulnerable to rupture at high stress rate (  ≥92.5EM Pa/s), resulting in TBI. The rupture of MTs results in accumulation of transported cargo in axonal swellings [7]. Fortunately, MT has the ability to recover to its original configuration through self-repair [13,68]. In addition, TBI may trigger chronic traumatic encephalopathy. Progressive axonal injury and structural degradation are considered as the classical features of chronic traumatic encephalopathy. These symptoms may be related to Parkinsonism and Alzheimer's disease [23,28,69].
From the perspective of biomedical engineering, this viscoelastic SLM is likely to become the core of the broader multi-scale model in time and space. As for the time, the computational window from milliseconds to years can explain axonal injuries in the gradual deterioration of the mechanism of time.
This can help identify the early markers of neurodegenerative disease and promote early treatment. As for the space, bridging from MTs to MAP tau can explain the molecular mechanism within the axon, which can provide new coping strategies to slow, block or reverse neurodegenerative disease [70]. This viscoelastic SLM may be a crucial step for shedding new light on the complicated interactions in MT-tau transient.
Though the viscoelastic SLM provides a new method to predict the axonal injuries under comparatively high stress rate over a typical timescale from microseconds to decades of seconds, the current study is not without limitations. A major limitation of this work is the choice of the Kelvin-Voigt model, which cannot represent stress relaxation. Surely if the axon is deformed and held at the final deformation, stress relaxation occurs in the tau protein. In order to accurately predict the underlying pathophysiology of DAI, future work will seek to improve the accuracy of the model by relaxing some assumptions and considering more of the related biomedical phenomena.

Conclusions
The aim of this study is to explore how the axon can simultaneously achieve a remarkable mechanical balance of high specific stiffness and toughness, and why the randomly staggered alignment microtubules are selected in the axon. Of course, viscosity is one of the significant physical attributes for axons [6,7], which will cause them to be sensitive to loading speed. In this study, we have extended the previous SLM by including Kelvin-Voigt viscoelastic behaviors for the tau protein of the axon subjected to transient loading. This Kelvin-Voigt viscoelastic SLM is an essential extension to the previous SLM which only considered the elastic or elasto-plastic behaviors under static or quasi-static conditions. Kelvin-Voigt viscoelastic SLM has been developed to elucidate the biomechanical behavior between microtubules and tau protein in the axon under transient loading. Theoretical closed-form solutions is presented to predict the stress transfer and failure forms while axon is stretched transiently.
According to the recent experiments [14,60] and results presented in this paper, the conclusions can be drawn as follows: (1) The maximum axial strain of tau protein is at both ends of the MT-tau interface length, and the minimum tau protein axial strain is at the center. The maximum and minimum axial strains of MTs are at the midpoint and both ends of the entire MT length, respectively.
(2) Axonal failure mechanism may be different under different loading conditions, such as microtubule rupture, detachment and detachment-rupture. In the process of MT-tau protein interface failure, tau protein serves as a molecular switch. Long microtubules are more vulnerable to rupture at high stress rate, yet short microtubules are likely to detach from the microtubule bundle under large deformations.
(3) Discontinuous, stiffest and randomly staggered alignment microtubules weakest-linked by tau proteins enable axon to withstand large deformations and have the advantage of dissipating energy. In this case, the forms of axonal mechanical failure is multi-level resulting in protecting the axon from overall damage. Therefore, composites with randomly staggered element can simultaneously achieve an outstanding mechanical balance of high strength and high toughness.