1 Introduction

Traumatic brain injury (TBI) has become a major public health problem, affecting more than 2.5 million individuals in the US^[1] and killing 15 per 100 thousand people in European countries each year^[2]. Numerous studies revealed that the mechanical breaking of axonal microtubules (MTs) in the white matter region leads to diffuse axonal injury (DAI), which is one of the most common pathologies of TBI. Indeed, DAI mainly results from overstretch or torsion of axon bundle through an external mechanical force such as rapid acceleration and deceleration, rapid rotation, contact forces, and penetration of a projectile^[3-4], causing a prolonged coma, patient's vegetative state, and subsequent death^[5]. Unfortunately, today's methods of imaging, such as conventional magnetic resonance imaging (MRI) and computed tomography (CT), are unable to detect the microscopic diffuse injuries to axons^[6-8]. Therefore, computational modelling is considered as an efficient tool to simulate the mechanical behaviour of brain tissue subject to external forces for substantially improving contemporary understanding of the underlying pathology and molecular mechanisms of axonal injury.

A single long slender axon and multiple shorter dendrites are specialized to transmit signals and receive signals, respectively^[9], as shown in Fig. 1(a). Histologically, axonal cytoskeleton is made up of MT bundles cross-linked by MT-associated protein tau^[9-10] (see Fig. 1(b)). It can be noted from the cross-sectional electron microscopy image (see Fig. 1(c)) that the MT bundle is generally hexagonally aligned^[11]. Parallel arrangement of MTs, acting as intracellular transport tracks and neuronal structural element, is the stiffest component of axon with the main structural function of reinforcing the neuronal axon^[12-15].

According to the recent experiments by means of atomic force microscopy, the shear modulus of MTs is 2-3 orders of magnitude smaller than Young's modulus^{[10, 17]}. Hence, the shearing effect may be the most critical factor leading to the underlying pathophysiology of DAI. The brain stem torques quite easily, which results in torsion of the nerve fibers in the core of the brain^[18]. In the view of animal rotational acceleration experiments^[19-21], rotational acceleration is suggested as the principal mechanical force causing DAI^[22-23]. Regardless static deformation or free vibration, the torsion mode is one of the principal shapes of failure modes of MTs^[24-26]. Practically, torsion of filaments is such a tricky experiment^[27-29] that computational modeling is an efficient method to study mechanical behavior of MT bundle and individual MTs under torsion. Buckling behavior of individual MTs under torsion was studied by Yi et al.^[30], and they found that a critical buckling torque is about 0.077 nN· nm. Chełminiak et al.^[31] developed a rigorous analysis of the MT elastic deformations in terms of the torsional degrees of freedom using the helix-based cylindrical structure of this biopolymer. A discrete bead-spring model was employed to analyze mechanical behavior of MT bundles, their relaxation behavior, and the influence of tau protein spacing under pure torsion^[3]. A finite element model was supposed to be an efficient tool to investigate the nonlocal torsional frequency and wave response of individual MTs^[32], or analyze the brain injury caused by rotational acceleration^{[6, 33]}.

Owing to its viscoelastic nature, axon usually exhibits strain or stress rate-dependent breaking. Axon bundle can withstand strains of over 100% and recover to its original configuration without evidence of damage^[34]. Nevertheless, failure of axons appears when rapidly stretched over a typical timescale of milliseconds^[35]. Thus, it is essential to construct a viscoelastic model to characterize the mechanical response of MT bundles under a dynamic loading. Generally, MTs are fundamental load-carrying element, and the tau protein is responsible for deformation. Hence, MTs and the tau protein are treated as elastic elements and viscoelastic springs, respectively.

Although some studies revealed the mechanical properties of MTs and the tau protein, the molecular mechanisms of axonal failure under rotational acceleration remain poorly understood. Noteworthily, chronic traumatic encephalopathy will continue to develop for decades after the initial injury; the more severe the original insult, the longer the survival time and the greater the severity of neurodegeneration. At the molecular scale of axon, random staggering of discontinuous MTs packed together by connective tau protein in a unidirectional nanobiocomposite. The aim of the present study is to develop an MT-tau dynamics model to explore the inherent mechanical vulnerability of axons in TBI, helping predict the underlying pathophysiology of DAI under the dynamic torsion condition. To investigate the micromechanical properties of axon both in space and time, an extended tension-shear chain model with Kelvin-Voigt viscoelastic tau proteins is constructed to account for the deformation and stress transfer between MTs and tau protein in the axon bundle under the dynamic torsion condition.

2 Model formulation and solutions

The tension-shear chain model was originally proposed by Gao et al.^[36] in an attempt to elucidate the path of load transfer in the mineral-protein biocomposites where the mineral platelets carry tensile load and the protein transfers load between mineral crystals via shear (see Fig. 2(a)). In the present study, we extend the tension-shear chain model to account for non-uniform staggering of axonal MTs in a unidirectional composite structure cross-linked by viscoelastic tau protein, where the MTs carry shear stress and the tau protein transfers load between neighboring MTs via tension under an applied torsional load, as illustrated in Fig. 2(b).

For a randomly staggered structure with an arbitrary overlap length between neighboring MTs, all MTs are assumed to be parallel aligned along the x-direction. They are hollow cylindrical structures with the outer radius R_O, the inner radius R_I, and Young's modulus E_M, which are assumed to be linear elastic materials. Neighboring MTs overlap nearly half of their length with an overlapping length donated by L.

It was reported that axon bundles have viscoelastic material behavior^{[1, 38-40]}. A Kelvin-Voigt viscoelastic model (see Fig. 3(c)) containing a spring with stiffness K and in parallel with a dashpot with the viscosity μ is applied for tau protein in an attempt to characterize the mechanical response of axon to dynamic torsional loading according to the recent experiments^[41],

(1)

where the tau protein dashpot timescale η =μ/K, μ is the viscous coefficient of tau protein, and the "overdot" above a symbol represents the derivative of the symbol with respect to time.

The torque T(x, t) leading to the elongation of tau protein can be converted to the normal stress σ_M(x, t) uniformly acting over the entire MT cross section (see Fig. 3(e)), which can be written as

(2)

where d_T is the center-to-center distance between the neighboring tau proteins.

For an MT bundle with hexagonal configuration (see Fig. 3(b)), if the rotations of the neighboring MTs are identical, there is no relative rotation between them, which implies that the tau protein has no deformation. If the rotations of the neighboring MTs are different from each other, a relative rotation occurs accompanied by the elongation of tau protein. Let us introduce the tau protein elongation δ(x, t), which is related to the difference of the individual rotation angle of neighboring MTs at the cross section x, as illustrated in Fig. 2(b). Hence, we can obtain

(3)

where γ is the angle between F₁ and F₂, and ϕ₁ and ϕ₂ are the rotation angles of the 1st and 2nd MTs, respectively. Additionally, note that the first and second derivatives with respect to x of δ(x) are denoted by δ' and δ", respectively. Other functions follow the same rule.

The circumferential relative rotation angle between the neighboring MTs will induce a shear stress acting over the MT cross section, as illustrated in Fig. 3(e). In the view of the classical torsion theory, the shear stress τ₁ and τ₂ over the 1st and 2nd MTs cross sections can be expressed as follows:

(4)

Obviously, the local torsional equilibrium of the neighboring MTs requires

(5)

where the torque applied to the 1st and 2nd MTs can be derived from T₁=F₁R₁ and T₂=F₂R₂, respectively. It is noted that T'₁ =-T'₂ as implied by the second term in Eq.(5). Therefore, the torque, which is caused by the circumferential shear stress, is identical in quantity.

We assume that the shear stress distributes uniformly over the cross-sectional area of the 1st and 2nd MTs. Thus, Eq.(5) can also be expressed as

(6)

where S=π(R_O²-R_I²) is the MT cross-sectional area.

Combining Eqs.(4) and(6) yields the equilibrium of the infinitely small section dx in the neighboring MTs as follows:

(7)

With the third term in Eq.(3) and the second term in Eq.(7), the following equations can be obtained:

(8)

Consider the equilibrium of an infinitely small section dx of the neighboring MTs in the circumferential direction of the bundle as illustrated in Figs. 3(a) and 3(d). We can obtain the 1st and 2nd MT equilibrium equations 2(R_O-R_I)σ_M(x, t)dx=Sdτ₁(x, t) and 2(R_O-R_I)σ_M(x, t)dx=-Sdτ₂(x, t), respectively, which can also be written as

(9)

Letting

(10)

and combining Eqs.(1)-(10), we can obtain the partial differential equation for the tau protein elongation δ,

(11)

The derivative of MT rotation angle with respect to x along the x-direction can be derived by combining the second term in Eq.(3) and the first term in Eq.(7),

(12)

with the boundary conditions ϕ'₁(0, t)=0 and ϕ'₁(L, t)=T cosγ/(GSR₁²). According to Eq.(12), the boundary conditions can be alternatively expressed as

(13)

while the initial condition can be expressed as follows:

(14)

where L_C is the characteristic length over which the stresses are transferred between the MT and the tau protein. The tau protein elongation can be derived by solving the nonhomogeneous partial differential equation(11) with the nonhomogeneous boundary conditions(13) and the initial condition(14),

(15)

It can be seen from Appendix A that, C₀ and C_n are determined by the initial condition of the tau protein elongation. If there is no elongation in the beginning, C₀=C_n=0. D₀ and D_n are derived from the force applied at the midpoint of MT. After the solution of tau protein elongation is obtained, the expressions for the shear stress(strain) distribution over the MT cross section and normal stress(strain) in the tau protein can be easily derived from the above equations.

The mechanical properties of axon with various viscosity coefficients or spring constants can be easily obtained according to the closed-form solution to this problem. In order to clarify the relationship between the applied stress rate and material properties of tau protein(the spring constant K and the viscosity coefficient μ), the time and spatial coordinate systems are rescaled as Y=t and X=x/L. In a similar manner, the variables in the governing equation(11) can be normalized as ∆(X, Y)=δ/L, , . In terms of the rescale variables, Eq.(11) can be expressed as

(16)

Equation(16) shows that the tau protein with high viscosity has similar stress distribution and deformation with high applied stress rate in the axon when all the other parameters remain unchanged. In contrast, the tau protein with high spring constant has opposite tendency.

3 Related experiments and parameters

According to the governing equation(11), the mechanical interaction behavior of MT bundle depends on two parameters, namely, (L_C/ζ)² and 1/ζ², of which the boundary conditions are independent. Actually, the values of these two parameters are determined by the material and geometrical parameters of the axon apart from the MT-tau interface length L. Note that the initial condition and boundary conditions are independent of these two parameters. Indeed, the boundary conditions depend on the MT-tau interface length L and the torque T(t) applied at the midpoint of MT. Therefore, the mechanical properties of axon are determined by the MT effective stiffness, the MT geometrical dimension, and the viscoelasticity of tau protein between MTs, which is in excellent agreement with the ultrastructural point of view proposed by Ahmadzadeh et al.^[12].

Experimental data from numerous studies are employed to assign values to the parameters used in the present theoretical solutions. Typically, an axon bundle comprises 10-100 MTs per cross section^[42]. MTs in the brain are approximately 1 μm-10 μm in length, which is hollow cylindrical tubes averagely formed by 13 tubulin protofilaments assembled in parallel, with outer diameter 25 nm and inner diameter 14 nm^[10]. MT bundles connected by the tau protein in the axon have a typical surface-to-surface spacing of ~ 20 nm and are hexagonally packed^{[11, 39, 43]}. The MT ultimate shear strength is defined as the shear strain value when the MT breaks. Atomic force microscopy(AFM) makes it possible for direct mechanical measurement at the micro-level. A study of the anisotropic mechanics of MTs has shown that the longitudinal αβ-αβ tubulin dimer bonds along protofilaments are a few orders of magnitude larger than the lateral bonds between adjacent protofilaments. Actually, the shear strain of MT is the result of two possible contributions: horizontal and vertical shear strains, which stem from distortion of inter- and intra-protofilament bonds, respectively. Intra-protofilament bonds are much stronger than inter-protofilament bonds, as reported in the literature^[10]. MT shear stiffness of torsion originates from the distortion of intra-protofilament bonds ranging from 1 kPa to 1 000 kPa^{[30, 44]}. In addition, buckling behavior of individual MT under torsion was studied by Yi et al.^[30], and they found that a critical buckling torque is about 0.077 nN· nm. This result was generally accepted by other researchers^{[3, 31]}, and the MT rupture shear strain can be calculated about 50%.

The other key parameters in our model are geometrical and material parameters of the tau protein. About 95% of tau filaments are in the form of paired helical filaments with a radius of 4 nm-10 nm^[45]. The tau protein center-to-center spacing ranges from 20 nm-40 nm according to a statistical analysis given by Hirokawa et al.^[46]. In accordance with the mean slopes of the force-extension curve, the spring constant of tau protein is estimated to be about 0.25 pN/nm, which is comparable to the stiffness of related proteins reported in the literature^[47]. Based on the experimental data^[41], Ahmadzadeh et al.^[12] estimated the viscoelastic parameter η = 0.35 s by using Bell's equation^[48]. In addition, interactions in hTau40 ruptured at peak forces 300 pN under quasi-static loading conditions observed by Wegmann et al.^[41]. Thus, tau protein breaking tensile strain is estimated in exceed of 40%.

After a comprehensive consideration from reference configuration, the geometrical and material parameters used in this study are obtained originally from experiments and estimations, which are listed in Table 1. For these parameters, we obtain the characteristic length from Eq.(10), L_C=32.84 nm, about 2 orders of magnitude smaller than the MT length.

4 Results and discussion 4.1 Process of axonal deformation

It can be noted from the dimensionless governing partial differential equation(16) that, the tau protein elongation δ(x, t) is determined by μ/K and L_C/L. The characteristic length, L_C = 32.84 nm, is much smaller than the MT length. According to the constitutive relation(1), the tensile force acting on the entire MT is not only closely relevant to the elongation of tau protein, but also influenced by the tau protein elongation rate. A unit cell, subject to the constant shear stress rate =20SR₁ kPa/s, is firstly studied. The 2nd MT is fixed at the end x=0, and the tensile stress is acting on the 1st MT at the other end x=L(refer to Fig. 3(d)). The distributions of tau protein elongation and MT shear strain at different time are shown in Fig. 4.

It can be noted from Fig. 4 that the values of the tau protein elongation δ and shear strain γ in MTs are all continuously increasing, while the applied shear stress at the end x=L is applied. The elongation of tau proteins located at both ends of the interface length are extremely large(see Eq.(15) and Fig. 4(a)). Generally, shear stress transfers from the loaded end towards the interior point of the MT-tau interface in the process of the dynamic torsion. At the loaded end, the MT shear strain is the maximum, and the minimum MT shear strain appears at the force-free end, which is equal to zero(see Fig. 4(b)).

According to the simple and elegant governing equation(see Eq.(11)) and the axonal deformation in the loading process(see Fig. 4), this MT-tau protein dynamics model is likely to become the core of the broader multi-scale computational model in time and space. As for time, the expansion of the simulation can explain several tau protein diseases in the gradual deterioration of the mechanism of time, which can help identify early markers of neurodegenerative disease and promote early diagnosis and treatment. As for space, bridging from the molecular level to the whole brain can explain how the macroscopic forces in TBI translate into cell and subcellular depletion, which can guide the design of new strategies to slow and block neurodegeneration^[50]. This computational model may be a significant step for shedding light on the complex micromechanical interactions between axonal MTs structure, tau protein, and mechanical dynamics.

4.2 Effect of MT-tau interface length on axonal deformation

Recently, we know that the interface length plays a significant role in mechanical behavior and stress transfer in fibrous structures^[51]. Two different MT-tau interface lengths subject to different loading rates are taken into account. The shorter one L = 328.4 nm, and the longer one L = 2 μm. The distributions of tau protein elongation and MT shear strain are plotted in Fig. 5. In the present study, the extended tension-shear chain model taking into account the viscoelasticity is developed in an attempt to elucidate the MT-tau interface length effects on axonal deformation under dynamic torsional loading.

Figure 5 shows that the MT shear strain distribution around the midpoint remains constant, if the MT-tau interface is long enough. As the interface length is much longer than the characteristic length(L≫L_C in Eq.(16)), e.g., L = 2 μm, the tau protein elongation in the region 0.1 < x/L < 0.9 is close to zero(see Fig. 5(b)). Furthermore, the results imply that the longer MT-tau interface is, the more severe stress concentrates at the loaded end and the more likely tau proteins will rupture.

4.3 Effect of applied loading rate on axonal deformation

Owing to its viscoelastic nature, the applied loading rate has a great effect on axonal deformation and stress transfer(refer to Eq.(16)). In Subsection 4.1, we have discussed the deformation of MT and tau protein when axon was twisted by different stresses with an identical stress rate. We are now to examine the mechanical response of axon under the same stress but at different stress rates, as shown in Fig. 5.

Figure 5 shows that the tau protein elongation at both ends(x=0 and L) is larger at a slower loading rate, but is nearly zero in the interior region, when the axon is twisted under the same stress but at different stress rates(see Figs. 5(a) and 5(b)), in excellent agreement with experimental observations^[52]. When the loading rate is fast enough, the MT shear strain distribution of the surrounding midpoint remains constant. Moreover, the results indicate that, with the increasing strain rate, the shear strain in the MT at the loaded end increases due to contributions from the tau protein viscoelastic effects. In a similar manner, we can note from Eq.(16) that the tau protein with high viscosity μ or low stiffness K leads to large elongation of tau protein in the MT at both ends.

4.4 MT and tau protein failure

The applied shear stress rate, the MT-tau interface length, the neighboring tau protein spacing, etc, have a vital impact on the mechanical failure of MTs, according to Eqs.(11) and(16). The relationship between these three parameters and axonal failure is shown in Fig. 6.

4.4.1 Effect of applied loading rate and interface length on axonal failure

With the increasing applied stress rate, MTs are more prone to rupture due to the contributions from the viscoelastic effects of tau protein. However, tau proteins are more likely to suffer mechanical failure under slower applied loading rate(refer to Figs. 5 and 6), which results in detachment of MT from the axonal MT bundle, which is consistent with the experimental observations^[3].

When the maximum shear strains in the MTs are equal in magnitude, the tau protein elongation in the longer MT-tau interface is larger than that in the short one. Conversely, the MT shear strain in the longer MT-tau interface is smaller than that in the short one, while the value of the maximum tau protein elongation is identical. It means that both the MT and tau protein are more likely to rupture in the long MT-tau interface. These phenomena are quite different from the mechanical behavior of fibrous structures subject to the dynamic tensile force^[51], for the reason that the axonal MT suffers from bending deflection and shear deformation in the torsional process, as observed by Lazarus et al.^[3].

4.4.2 Effect of tau protein spacing on axonal failure

According to Eq.(10), the characteristic length L_C is dependent on the spacing between tau proteins, and also L_C increases with increasing d_T. By doubling the spacing of the tau proteins when compared with the previous subsection, the value of L_C increases to 46.44 nm. In this case, the shear strain in MTs as a function of tau protein elongation along the short and long MTs is plotted in Fig. 6.

Comparing with different spacings between tau proteins illustrated in Fig. 6, we can find that upon increasing the spacing between tau proteins, the threshold tau protein elongation for MT rupture generally decreases in all the cases.

In the present study, the tau proteins are assumed to uniformly distribute along the MTs. Actually, the distribution of tau proteins along the MT is random, to which the extended tension-shear chain model is equally applicable. The entire MT-tau interface length should be divided into several short segments according to the real distribution of tau protein along the MT, in such a way that deformation of MTs should satisfy the continuity conditions between adjacent segments. Therefore, the deformation and stress distributions of MT and tau protein are multi-level.

TBI, a non-degenerative and non-congenital insult to the brain caused by sudden external mechanical forces, is one of the most common public health issues. The rupture of MTs as the tracks for chemical cargo transport leads to accumulation of transported materials in axonal swellings^[33]. Fortunately, MTs have the self-healing properties. Therefore, MTs recover their mechanical stiffness through self-repair^{[15, 53]}. Moreover, the breaking of tau protein will accumulate in the axon, which may trigger progressive neurological degeneration^[54]. It is striking that these symptoms appear to be associated with many other neurodegenerative diseases, including Parkinsonism and Alzheimer's disease^[55]. Thus, the understanding of the mechanism of the interaction between the axonal MT and tau protein is extremely important.

5 Conclusions

It is extremely difficult to experimentally measure the micromechanical properties of MT bundles under torsional loads^[27-28]. In this study, an extended viscoelastic tension-shear chain model is developed to computationally investigate the mechanical behavior of MT bundles under external torsions. To the best of the authors' knowledge, there are few pathological reports available for characterizing the mechanism of axonal MT bundles resisting torsional force, yet rotational acceleration will surely cause DAI and neurodegenerative diseases. Obviously, a mere computational model, without validation by experimental or other computational studies, could not provide comprehensive and accurate results on this subject. However, the accuracy of the MT bundle model was assessed in two studies on tension of MT bundles^{[12, 34]}, where numerous experimental results were used for validation. This study will lay the foundation for future computational and experimental studies on axonal cytoskeleton to predict the damage in the brain and understand the underlying pathophysiology of DAI. In summary, conclusions can be drawn as follows:

(ⅰ) The large elongation of tau protein appears at both ends. The MT shear strain is the maximum at the loaded end, and the minimum MT shear strain appears at the force-free end.

(ⅱ) The longer MT-tau interface is, the more severe MT shear stress concentrates at the loaded end and the tau protein normal stress concentrates at both ends. Hence, both the MT and tau protein are more likely to rupture in the long MT-tau interface.

(ⅲ) With the increasing strain rate, the shear strain in the MT at the loaded end increases due to the contributions from the tau protein viscoelastic effects. MTs are more prone to rupture at a fast loading rate(or the viscosity is very high). However, the tau protein is more likely to suffer mechanical failure at a slow loading rate(or the viscosity is small).

(ⅳ) The threshold tau protein elongation for MT rupture generally decreases with the increase in the spacing between tau proteins.

Particularly, while the emphasis of the present study is on the dynamic micromechanical behavior of axon, this viscoelastic tension-shear chain model is equally applicable to other soft tissues(ligament, muscle, tendon, etc.) possessing the similar hierarchical microstructure. Furthermore, the present theoretical method can also be used to explain the deformation and stress transfer in the hard tissues(bone, tooth, etc.) with approximately zero viscosity in the governing equation(11).

Appendix A

The tau protein elongation δ(x,t) is obtained by the theoretical derivation via the following steps.

(ⅰ) The solution δ(x,t) is split into two parts V(x,t) and W(x,t),

(A1)

W(x, t) carries the burden of the nonhomogeneous boundary conditions. Therefore, W(x, t) is assumed to satisfy the boundary conditions,

(A2)

The simplest function used for W(x, t) is

(A3)

(ⅱ) Substituting Eqs.(A1) and(A3) into Eqs.(9)-(12), we can obtain the nonhomogeneous partial differential equation with homogeneous boundary conditions and nonhomogeneous initial conditions, which can be expressed as

(A4)

where

(A5)

(ⅲ) Equation(A4) can be divided into the homogeneous partial differential equation V₁(x,t) with nonhomogeneous initial condition and the nonhomogeneous partial differential equation V₂(x,t) with homogeneous initial condition,

(A6)

(A7)

(ⅳ) The theoretical solution to the problem(A6) can be derived by using the method of separation of variables,

(A8)

where

(A9)

(A10)

(ⅴ) We propose to solve this problem(A10) by using the method of eigenfunction expansion via expanding the unknown solution V₂(x, t) and nonhomogeneous term of partial differential equation H(x, t) in a series of the related homogeneous eigenfunctions,

(A11)

where