Annaclaudia Montanino
Marzieh Saeedimasine
Alessandra Villa
Svein Kleiven- 1Division of Neuronic Engineering, Royal Institute of Technology (KTH), Stockholm, Sweden
- 2Department of Biosciences and Nutrition, Karolinska Institutet (KI), Stockholm, Sweden
Traumatic brain injuries are a leading cause of morbidity and mortality worldwide. With almost 50% of traumatic brain injuries being related to axonal damage, understanding the nature of cellular level impairment is crucial. Experimental observations have so far led to the formulation of conflicting theories regarding the cellular primary injury mechanism. Disruption of the axolemma, or alternatively cytoskeletal damage has been suggested mainly as injury trigger. However, mechanoporation thresholds of generic membranes seem not to overlap with the axonal injury deformation range and microtubules appear too stiff and too weakly connected to undergo mechanical breaking. Here, we aim to shed a light on the mechanism of primary axonal injury, bridging finite element and molecular dynamics simulations. Despite the necessary level of approximation, our models can accurately describe the mechanical behavior of the unmyelinated axon and its membrane. More importantly, they give access to quantities that would be inaccessible with an experimental approach. We show that in a typical injury scenario, the axonal cortex sustains deformations large enough to entail pore formation in the adjoining lipid bilayer. The observed axonal deformation of 10–12% agree well with the thresholds proposed in the literature for axonal injury and, above all, allow us to provide quantitative evidences that do not exclude pore formation in the membrane as a result of trauma. Our findings bring to an increased knowledge of axonal injury mechanism that will have positive implications for the prevention and treatment of brain injuries.
Introduction
Traumatic brain injury (TBI) is defined as “an alteration in brain function, or other evidence of brain pathology, caused by an external force” (1). In 2013, ~2.8 million TBI-related emergency department visits, hospitalization, and deaths occurred in the United States (2). In a recent study reporting the epidemiology of TBI in Europe, Peeters et al. analyzed data from 28 studies on 16 European countries and reported an average mortality rate of ≈11 per 100,000 population over an incidence rate of 262 per 100,000 population per year (3). Diffuse axonal injury (DAI), a multifocal damage to white matter axons, is the most common consequence of TBIs of all severities including mild TBIs or concussions (4). Invisible to conventional brain imaging, DAI can only be histologically diagnosed and its hallmark is the presence of axonal swellings or retraction balls observable under microscopic examination (5).
Considerable research effort has been put into the understanding of the primary effects of trauma onto the neurons. To define an axonal injury trigger, nervous tissues and neuronal cultures have been subjected to dynamic loads leading to several hypotheses regarding the cell injury mechanism. Mechanoporation (i.e., the generation of membrane pores due to mechanical deformation) of the axolemma has been historically put forward as primary axonal injury mechanism (6–11). Although such a mechanism was initially proposed to affect the somatic plasmalemma, it was more recently associated with both the somata and neural processes (12, 13). Assessing the occurrence of such a mechanism is particularly important in cases of mild/moderate TBIs, where cells' abnormalities—but not immediate death- are coupled with functional alterations (14). Disruption of the axolemma has been inferred observing changes in membrane permeability due to mechanical stretch supposedly leading to increased intra axonal calcium concentration, however this has been questioned by some studies (15–17).
Some experimental studies using oriented cell cultures have instead put forward microtubule disruption as cell injury trigger having observed cytoskeletal damage and tau protein accumulation at the site of axonal swelling (18, 19). These studies have further motivated the development of models (both analytical and computational) having the microtubule bundle as a focus (20–25). In our previous work, however, we have shown with a finite element model of the entire axon that, especially when including detachment of tau protein elements, deformations in the microtubules did not pass the suggested microtubule failure thresholds (26). On the contrary, we found that, as a result of axonal stretching and of microtubules distancing, very high strain localization formed on the axonal membrane. However, we were unable to say whether these could lead to axolemmal damage or not.
Up to date, studies assessing poration thresholds—both experimental and computational ones—have not been axon-specific and therefore might not reflect axolemmal behavior. Previous experimental studies have in fact focused on areal deformations of red blood cells (RBCs) (27). Molecular simulations studies have investigated mechanoporation by deforming two-to-four component lipid bilayer biaxially (28–30). Given the slenderness of the axons, however, during stretch-injury one can expect the axolemma to deform mainly along the axial direction. Not only loading mode specificity, but also material specificity might also play a key role in the understanding of axonal injury. Among others, lipid composition is known to affect membranes' biomechanical properties (31). For example, experimental studies have established areal failure thresholds (εA < 10%) based on RBC's or giant lipid vesicles' membranes (32, 33). It has however been shown that these thresholds do not apply to the axons' membrane: when undergoing osmotic shocks—which induce biaxial deformation of the axolemma- in fact the axolemma could sustain areal strains εA > 20%) without rupturing (34).
Another peculiarity of the axonal membrane is that both in unmyelinated and myelinated axons, protein channels span the lipid bilayer, making possible the propagation of the action potential. In particular, voltage-gated sodium channels are the main responsible for the continuous and saltatory conduction in unmyelinated and myelinated axons, respectively. The density of sodium channels in the axonal membrane has been reported to be in the range 5–3,000 channels/μm2. Lower densities of sodium channels are found in unmyelinated axons, while higher densities are observed in nodal portions of myelinated axons (35, 36). These channels have also been proposed to influence the injury response (37). Rigid inclusions embedded in a lipid bilayer model have been shown to facilitate membrane disruptions while stiffening the membrane-inclusion system (38). In this study, unmyelinated axons are considered. These axons are, in fact, not only as numerous as myelinated ones in the human brain (19), but have also been shown to be more vulnerable to mechanical strain than their myelinated counterpart (39, 40). These motivations have so far justified the usage of unmyelinated axons as experimental injury models. Modeling an unmyelinated axon allows us to gain a better insight of the basic axon injury mechanism, irrespective of more complex and ill-defined boundary conditions—such as those that could be induced especially at the paranodal junction by the attachment of myelin or other surrounding tissue.
It is apparent that, in order to describe the axonal injury mechanism, both axon-specific boundary conditions and material properties should be taken into account. Therefore, to provide mechanical insights into the initiation of axonal damage, we combined a finite element (FE) model of a generic distal portion of an unmyelinated axon and a molecular-based lipid bilayer model with fixed number of particles. In particular, the axonal model was utilized to simulate typical stretch injury scenarios (Figures 1A,D). As a result of cellular-level deformation, local deformations happening at the cortex level (Figures 1B,E) could be extracted. These local deformations were then used as input for molecular simulations of the axonal lipid bilayer (Figures 1C,F). In this way, membrane permeability or poration could be quantified in dependence of applied axonal strain and strain rates.
Figure 1. Schematic representation of the combined modeling approach. (A) Finite element axonal model coated with a slightly transparent layer symbolizing the lipid bilayer (virtual). A quarter of the model was removed to reveal the arrangement of the inner structures. A microtubules bundle (solid blue hollow cylinders cross-linked by beam-like tau proteins) is located centrally, while the rest of the axonal space is filled with a dense beam-meshed neurofilaments network (yellow). The cytoskeleton is wrapped in a thick shell layer representing the axonal cortex (gray). (B) Single undeformed cortex element coated with a virtual lipid bilayer linked to the corresponding undeformed molecular coarse grained model of the lipid bilayer. (C) The lipid content resembles the plasma membrane: (PC lipids are in blue, PE in yellow, CHOL in green, SM in orange, PS in purple, GM in red, anionic lipids in pink, and the rest of lipids in gray, water is not visualized for clarity). (D) Deformed axon finite element model (axonal deformation, εaxon). (E) Deformed cortex element (localized deformation: εx = 0.34, εy= 0). (F) Visualization of the induced pore upon lipid bilayer deformation.
Methods
Finite Element Framework
FE head models are commonly used by the scientific community, despite their level of idealization, for the explanation they can provide of traumatic phenomena. Through FE simulations it is in fact possible to simulate traumas and extract quantities that would otherwise be inaccessible via an experimental setup. Similarly, the FE method has been used to capture the mechanical behavior of cells and subcellular components (41–47). More recently, FE models of the axon have been proposed to shed a light on the mechanism behind axonal injury (25, 26, 48–50).
In the present study, a previously published and validated axon FE model was used (Figure 1A), which is a 8 μm long representative volume (RV) of an axon consisting of three main compartments: a microtubule (MT) bundle, the neurofilament (NF) network and, finally, the axolemma-cortex complex wrapping the entire structure (26). The axon model diameter (1.15 μm) was chosen to fall in the range of those from experiments against which the model was originally validated against (51–53). The choice of its length was instead dictated by considerations related to the average MT length. Given the average continuous length of MTs (4.02 μm) (54), to contain a unit length of the microtubule bundle the length of the axon RV was set to 8 μm. Given its cross-section and the proposed MTs densities, 19 rows of MTs were included in the axon RV, each containing two randomly placed discontinuities (55–57).
In the MT bundle, solid elastic MTs are cross-linked by viscoelastic tau proteins, which are modeled as discrete beams. These are assigned a failure threshold so that they stop carrying the load as soon as they reach double their length (20). The NF network is modeled as a dense mesh of viscoelastic beams filling up the axonal space and anchoring the MTs to the plasma membrane (58). A detailed description of our literature-supported modeling choices can be found in the Supplementary Material and in our previous publication (26).
In a study focused on the assessment of RBC membrane mechanical behavior, Evans et al. defined the lipid bilayer as “along for the ride” when the cell is deformed, meaning that the cortex (or “matrix”) is the one resisting the applied deformation and influencing the structural response (59). Assuming this to be true also for the axon, in the current study, the axonal wrapping layer was assigned properties in line with those of the sole axonal cortex. The axonal cortex has a unique periodic organization made of rigid actin rings connected by flexible cylinders of spectrin (60). In our model the cortex is represented by a layer of fully integrated shell elements with a thickness of 50 nm (61). These elements were assigned a shear modulus G = 0.0016 MPa and a viscosity η = 105 Pas (62, 63). Mechanical properties of all the model components can be found in Table 1. In a recent study, fluid friction was observed at the membrane-cortex interface with friction forces being proportional to the relative speed of the two layers (61). Considering the axonal loading mode chosen for this investigation (i.e., pure stretching), however, the relative movement between the two layers can be assumed to be minimal. Hence, we assumed the lipid bilayer as fully tied to the cortex: any deformation happening in the cortex plane would be directly transferred to the lipid bilayer. Therefore, in the current study an explicit FE representation of the membrane was avoided (hence the “virtual” representation in Figure 1) and deputed instead to a molecular dynamics description. In practice, cortical deformations resulting from axonal injury simulations were extracted and fed into coarse grained (CG) membrane or membrane-protein systems.
Table 1. Axon model material parameters.
Axonal Injury Simulations
To study axonal injury at the single-cell level it was deemed appropriate to consider the axonal behavior under uniaxial deformation. It is generally accepted that strain is the main mechanism behind axonal injury (66). Although at the tissue (i.e., brain) level strain presents itself as a three dimensional (3D) tensor, to enforce experimental control, it is common to apply only one deformation at a time, this being of compressive, tensile or shearing nature. It has previously been shown with neuronal cultures that different loading modes might actually lead to different mechanisms of injury (67). Specifically, in that study, Geddes-Klein et al. found that stretching a primary cortical neuron culture uniaxially or biaxially (simultaneously in two perpendicular directions) yielded different mechanisms behind the influx of calcium. However, it should be considered that when deforming a network of cells—which are not embedded in a matrix-, the network's almost 1D structures, the single axons, will first reorientate in the direction of the load and then mostly stretch.
In an effort to investigate a general axonal behavior, rather than a particular one, 10 different FE axon models were produced by altering the MT bundle geometry of the baseline model. More specifically, every model was generated by randomly moving the MTs discontinuities locations, while keeping the average MTs length of the original model. Symmetry conditions were enforced on a quarter of a model and simulations were run with the latter condition.
Each axon FE model was subjected to a displacement-controlled uniaxial tensile deformation (Figure 1D). Namely, each model was stretched up to a global strain εaxon = 30%. This deformation range was chosen to cover all the axonal injury thresholds proposed so far in the literature at either cellular or tissue level (68–72). The deformations were applied at strain rates 1, 10, 20, 40/s, which are characteristic of axonal injury experiments (7, 9, 18). It is important to note that our axon FE model represents an unmyelinated axon and hence can be compared against in vitro injury model, which typically do not include myelinated axons. All simulations were performed in LSDYNA using an implicit dynamic solver. The implicit scheme was a necessary choice dictated by the element dimension, while the dynamic regime was chosen to be able to capture strain-rate related inertia and material effects. Subsequently, for each simulation of 1, 10, 20, and 40/s strain rates the 1st and 2nd principal strains (εx, εy) in the cortex plane were extracted as function of axonal strain.
Setting Up the Membrane Molecular Model
The membrane was modeled as a molecular-based lipid bilayer (two leaflets of lipids) and was described by a CG model where groups of atoms (3–4 heavy atoms) are united into beads which interact with each other by means of empirical potentials. The MARTINI2.2 force field together with the non-polar water model was used (73–75). The axolemma is characterized by having phosphatidylcholine (PC), phosphatidylethanolamine (PE), cholesterol (CHOL), glycolipid, [and for some studies also sphingomyelin (SM)] as the main components (Table S1). To mimic the axolemma's lipid composition, the mammalian plasma membrane model deposited on MARTINI webpage (http://www.cgmartini.nl/) was used (76). In the membrane model, 63 different types of lipids were distributed asymmetrically between two leaflets (Figure 1C). The outer leaflet has a higher level of saturation of the tails and contains PC (36%), PE (6%), CHOL (31%), SM (19%), glycolipid (GM) (6%), and other lipids (2%). The inner leaflet, which has a higher level of polyunsaturation, contains PC (17%), PE (25%), CHOL (29%), SM (9%), phosphatidylserine (PS) (11%), anionic lipids phosphatidylinositol (PIP) (2%), and other lipids (7%). Note that the reported values are in mol%. Full details of the membrane model can be found in the work by Ingólfsson and coworkers (76). The bilayer (containing a total of 6,700 lipids) was placed in a cubic box (42*42*12 nm) and solvated by about 100,000 CG water beads and NaCl was added to mimic the ionic strength at physiological condition (150 mM NaCl).
Setting Up of Membrane-Protein Molecular Model
Sodium channels are abundant in the axonal membrane (77, 78). When the present study was conducted, enough information was available to model solely sodium channel protein type subunit alpha (Nav1.1) (79). Although, in fact, the 3D structure for the homo sapiens Nav1.1 has not been experimentally resolved, its encoding gene is known (SCN1A gene) (80). Thus we built the 3D structure using homology modeling and Phyre2 server was used (81). Among the available sodium channels, the best sequence alignment (100% confidence and 48% sequence identity) was found with putative sodium channel from American cockroach, NavPaS (PDB ID: 5X0M) (82). The cryo-electron microscopy structure of NavPaS has been used as template. The missing N-terminal and C-terminal domains, 1–2 and 3–4 linkers were modeled using the software Modeler (83).
The so obtained 3D structure for Nav1.1 is composed of a single polypeptide chain that folds into four homologous repeats (Figure 2), each one containing six transmembrane segments. The protein was oriented in the lipid bilayer using Positioning of Proteins in Membrane (PPM) method (84). The structure was first minimized in the PC bilayer at the atomistic level [using CHARMM36 force field (85)] with 5,000 steps of steepest descent method and then equilibrated for 200 ns.
Figure 2. (A) Three dimensional structure of Nav1.1 protein embedded in the lipid bilayer. Protein is shown by secondary structure elements (α-helix in magenta, β-sheet in yellow, turn in cyan, and coil in black). Phosphate atoms of lipid head group are shown in blue. (B) CG model of Nav1.1 protein in the plasma membrane model. Protein is shown in black VDW representation. For lipid color code, see Figure 1.
The obtained atomistic structure was used to build the protein CG model. Martinize protocol together with ElNeDyn were used to convert the atomistic Nav1.1 structure into a coarse-grained model (86, 87). One protein was embedded in the 42 × 42 nm2 CG membrane (Figure 2), that corresponds to a density of 567 ion channels/μm2.
Molecular Dynamics (MD) Simulations
MD simulations are widely used to provide a 3D description of lipid bilayers and proteins at the molecular level and their evolution in time, from which time-dependent and independent properties can be evaluated (88). All MD simulations were performed using the GROMACS simulation package, version 2016 (89). Simulations have been performed in NPT ensemble (constant number of particles, constant pressure, and constant temperature) and NPzAT ensemble (constant number of particles, constant pressure along z, Pz, constant bilayer area, and constant temperature). The temperature of the systems was kept at 310 K using velocity rescale thermostat with a time constant of 1 ps (90). Semi-isotropic Parrinello-Rahman barostat was applied to couple the pressure to 1 bar with a time constant of 12 ps (compressibility of 3e-4 bar−1) (91). Periodic boundary conditions were applied. A time step of 20 fs was used. The Verlet cutoff scheme was used. Non-bonded interactions were calculated using a cut-off of 1.1 nm (92). The reaction field potential was used to treat long range electrostatic interactions using a switching distance of 1.1 nm (93).
NPT simulations were performed to equilibrate the membrane structure in absence of deformation. To generate the starting structures for deformed lipid bilayers, we deformed the bilayer using unsteady stretching algorithm with a stretching speed c = 5 m/s. The so obtained deformed system was then simulated at constant bilayer area for 2 μs to generate an ensemble of configurations that describes the membrane at the selected strain. The simulations were extended up to 5 μs for strains larger than 30%. During the simulations an elastic bond force constant of 500 kJ·mol−1·nm−2 was applied to the protein atom pairs within a 0.9 nm cut-off (ElNeDyn) (87). That means that the protein is simulated as a semi-rigid body. The analysis was performed on the last 1 μs of production runs. The errors were calculated by block averaging over five blocks. All the figures were rendered using Visual Molecular Dynamics software (94). As criteria to detect pore formation in the lipid bilayer, we use the jump in the surface tension. After the detection of the pore, simulations were extended of 5–9 μs to check whether spondaneous sealing occurred in the simulation time scale (μs).
Permeability Coefficient
The transport of substances across a lipid membrane is a biological process of vital importance. Small molecules, such as water molecules or drugs, can be transported into the membrane in a passive or active way. A passive transport proceeds via an entropy-driven, non-specific diffusion process of the molecule across the membrane. Permeability or leakage of a small molecule should give a measure of the structural stability of the membrane since it should reflect lipid packing at the membrane core. Membrane permeation can be described by a so called “solubility-diffusion mechanism” [see review by Shinoda (95) for details] and quantified by the membrane permeability coefficient, Pm. The Pm of a water molecule is proportional to its oil/water partition coefficient (Koil/water) (96, 97).
First, we calculated the membrane/water partition coefficient for water at equilibrium and at different strains. We do that by calculating the average number of water molecules (1 CG water bead = 4 water molecules) in the bilayers within a distance of 0.52 nm from the lipid tails on 1 μs molecular simulations. The partition coefficient for water is , where [solute] denotes the concentration of water molecules in the lipid bilayers and water solution, respectively. [water]mem is obtained by dividing the average number of water molecules by the lipid bilayers volume, while for the water concentration in water the experimental value (55.5 mol/L) was used. The membrane volume was corrected to account for the presence of the protein. Knowing that logPm is proportional to logKmem/water, we then estimated how membrane strain influences its permeability.
Results
In this study, axonal injury was simulated by deforming a FE model of the axons of a quantity εaxon. Cortex principal strains (εx, εy) were then extracted and applied in cascade to a molecular model of the plasma membrane and the structural outcome on the latter was observed. As evidenced by Figure 3—which shows the distribution of 1st Principal Green-Lagrange strains along the cortex shell-layer for one of the 10 tested models for 5% < εaxon <15%—when the axon is stretched (0 < εaxon <30 %), due to its composite nature, the deformation pattern along the axonal cortex is not homogeneous. Areas of strain concentration on the cortex evidently form throughout axonal deformation. The magnitude of this strain concentration is dependent both on the applied strain (different columns in Figure 3) and strain rate (different rows in Figure 3). It can be observed that, for example, an applied axonal strain as low as 5% corresponds to a local maximum strain which is twice as high. The more εaxon increases, the more this localization of strains in the cortex is pronounced. Moreover, on each column it is possible to observe, for the same axonal strain, the influence of strain rate (1, 10, and 40/s) on the local deformation level. Here it is possible to notice that higher strain rates do not lead to higher local maximum strains for εaxon < 10%. However, at higher axonal strains the intuitive order is reestablished: higher strain rates lead ultimately to higher cortex-level deformations (Figure S2).
Figure 3. Fringe plots showing the 1st Principal Green-Lagrange strain along the axonal cortex as a result of 5, 10, and 15% axonal strains in the first, second and third column, respectively. In the first, second and third row results for strain rates of 1, 10, and 40/s are reported. The range was set between 0 and 40% for visualization purposes.
To better understand this behavior, the average results of 10 different models are reported in Figure 4 together with the standard deviations for each tested strain rate. In this plot the cortex maximum 1st Principal Green-Lagrange strains (εx) are reported as a function of the applied axonal strains. What can be noted is, first, the increased standard deviation with higher rates. In addition, one can notice that for strain rates of 1 and 10/s the curves follow the same path until εaxon ≈ 7% and then diverge. Interestingly, as previously noted, the response for a strain rate of 40/s is on average lower than that of 1 and 10/s until εaxon ≈ 12% and εaxon ≈ 20%, respectively. Second Principal Green-Lagrange strains (εy) in the cortex plane were consistently found to be at least 6 orders of magnitude inferior to εx, therefore they were set to zero when input into MD simulations.
Figure 4. Maximum 1st Principal Green-Lagrange strains in the axonal corex as a function of axonal strain. Lines represent the average over 10 models and the shaded areas represent the standard deviations. Results are shown for strain rates 1,10, 20, and 40/s. In the second row (center), the means for the different strain rates are shown in the same plot to ease the visual comparison. Finally, in the second row (right), iso-strain curves are reported, which show the maximum 1st Principal Green-Lagrange strains in the cortex as a function of the applied strain rate.
To understand the effect of axonal deformation on the membrane integrity, we used MD simulations and a molecular based-lipid bilayer to describe the axonal membrane. First, to verify that the membrane models resemble the mechanical feature of natural membrane, the area compressibility modulus was estimated for small areal strain values (<0.05). A value of 303 ± 6 and 352 ± 3 mN/m, respectively in absence and presence of Nav1.1 protein, was obtained (Figure S3)—which is in agreement with the area compressibility modulus measured for RBC membranes (375 ± 60 mN/m) at 37°C by Waugh and Evan (98). The agreement allows us to move forward and use the model to investigate the effect of local deformation on the bilayer structure. A selected number of cortex strains εx (Figure 4) were input into the membrane and membrane-protein system, while the membrane model was kept undeformed in the y-direction (εy = 0).
We monitored the bilayer structural features (i.e., pore formation, interdigitation, and water permeability) at different local strains. Figure 5 depicts how the occurrence of pore formation is detected. Namely, this event corresponds to a jump in the surface tension. Figure 6 summarizes which molecular-level event is to be expected at each axonal strain εaxon (and related maximum cortex strain, εx) for a plasma bilayer model with and without embedded protein. The lipid bilayer can withstand the applied strain value of εx = 34%, corresponding to εaxon = 10–12%. At higher strains pore formation is observed. When Nav1.1 is embedded in the membrane, poration occurs at a larger cortex strain (εx = 47%), corresponding to εaxon = 12–15%.
Figure 5. Surface tension and corresponding lipid bilayer structure at εx = 0.34 as a function of simulation steps.
Figure 6. Axonal strains and corresponding maximum local strains (mean values over a family of 10 different FE axonal models) together with strain rates. *indicates the local deformations for which molecular dynamics simulations have been performed. Top panel (A) refers to the membrane model and bottom panel (B) to the membrane-protein model. ΔPm indicates how water permeability increases with respect to the equilibrium value (at local strain = 0). max corresponds to 40% in the membrane model and 51% in the membrane-protein system.
Figure 7 shows an example of the observed bilayer disruption in absence and presence of Nav1.1 protein (Movies S1, S2). We observed that the membrane rupture occurs at slightly larger local strain in presence of the protein. The presence of the Nav1.1 protein therefore not only renders the membrane more resistant to deformation but also to rupture. In addition, it was observed that before poration, the bilayers start to form interdigitated states (Figure S4). Interdigitation occurs at εx > 27% in absence of protein and εx > 34% in presence of protein, corresponding to εaxon =9–11% and εaxon = 10–12%, respectively (Figure 6). Interestingly, poration both in presence and absence of protein takes place in bilayer regions lacking ganglioside lipids (Figure 7). Ganglioside-type lipids are mainly found in the outer leaflet of the membrane and are known to make the lipid bilayer more resistant to deformation due to the interactions between sugar headgroups (30).
Figure 7. Lipid bilayers disruption. Left panels: snapshot at 680 ns at εx= 0.34 in absence of Nav1.1. Right panels: snapshot at 27 ns for εx = 0.47 in presence of Nav1.1. (in black). Top panels show the upper leaflets, low panels show the lower leaflets. For lipids color code, see Figure 1.
To quantify that the molecular models indeed reproduce membrane partitioning, we calculated the membrane/water partition coefficient for water at equilibrium (logKmem/water = −1.73). The presence of the protein has almost no effect on the partition coefficient (logKmem/water = −1.75). At the increase of εaxon, water partitioning between aqueous solution and lipid bilayers increases proportionally to εx (Figure S5). Given the aforementioned proportionality between the partition coefficient Kmem/water and the permeability Pm, the latter can be quantified as a function of strain. In particular, starting from equilibrium state, water permeability increases up to 40% (51% in presence of protein) before pore formation. Figure 6 coloring describes the change in water permeability for the the lipid bilayers with and without protein channel at the increase of εaxon for different strain rates.
Discussion
The biomechanics of TBI has been the focus of considerable research effort for many years now. Nevertheless, the underlying injury mechanism leading to cellular impairment is yet to be explained. The use of in-silico models can provide better understanding of the mechanical cues determining cells' condition immediately after a mechanical insult. In this study, for the first time, we were able to provide direct mechanical evidence of the possible onset of mechanoporation as a result of the mechanical insult. More specifically, bridging finite element and molecular dynamics simulations, we could determine at which level of axonal strain the lipid bilayer sustains local deformations high enough to induce pore formation.
Several studies so far have stressed the need of taking into account axonal direction and kinematics when studying injury at the tissue level both in experimental set-ups and computational models (99–101). It has indeed been assessed from tissue or cell cultures that deformations applied to a tissue do not necessarily coincide with those sustained by the embedded axons. In fact, in both tissue models -such as the optic nerve or organotypic slices- and cell cultures non-affine deformations might arise. Similarly, at the cellular scale, which is the focus of the current study, the applied axonal strains εaxon do not directly transfer to the axonal membrane, due to the axonal composite nature. Despite being often assumed as a homogeneous viscoelastic entity, the axon presents dramatic heterogeneities in deformation along its length when stretched (102). Evidence of this is also the localized axonal length increase (up to 65%) that was observed as a result of axonal stretch (19). The morphological response to stretch injury—such as the appearing of undulations and periodical swellings—clearly hints at a heterogeneity in the axonal response to tensile injury loads. This is confirmed by the localization of membrane strains that we observed with our model (Figure 3). This localization is directly related to the point of weakest connectivity of the microtubule bundle (24, 103).
Especially when considering that our axonal finite element model is just an 8 μm-long representative volume of the axon, one should imagine such a strain localization taking place periodically along a periodical repetition of our model. The periodical swelling, which is characteristic of axonal injury, could therefore potentially be explained by the heterogeneous deformations that are revealed by our model.
Recently, an experimental work has tried to quantify axonal membrane strains induced by the deformation applied to a cell culture as a whole (104). Differences between applied strains and membrane strains were reported. However, it must be noted that these discrepancies may mostly be due to the non-alignment between applied load and axons' directions within the culture. Moreover, membrane strains were calculated from the level of separation of fiducial membrane markers whose initial distance is not reported. Hence, a direct comparison between our results and these data cannot be made.
In our results maximum cortex strains clearly appear as a non-linear function of axonal strain (Figure 3). Less intuitive is however the relation with strain rate. Due to the viscoelastic (and hence rate-dependent) properties of tau proteins and of the cortex itself, at axonal strains lower than 12%, the cortex deforms less at strain rate 40/s rather than at 1 or 10/s. Similar observations were reported in previous studies based on numerical and analytical models of the sole microtubule bundle (24, 25, 103). In those studies, simulations of the microtubule bundle undergoing stretch revealed a higher resilience of the microtubule bundle at higher stretch rates. Whether this emerging property is reflected by experimental models of axonal injury has, to the best of our knowledge, not been established yet.
To assess the influence of the chosen material properties on maximum cortex strains and on this rate-related phenomenon, a sensitivity study was conducted with one of the ten geometries (Figure S1). The results indicate that, while neurofilament and cortex properties minimally affect cortex strains, tau protein viscosity does, particularly at rate 40/s. Changes in tau proteins' viscosity yield similar results at rate 1 and 10/s. However, at rate 40/s increasing or decreasing the viscosity of 100%, delays or anticipates, respectively, the crossing of the mechanoporation threshold. A purely elastic behavior for these elements (null viscosity) reestablishes the intuitively expected order: cortex maximum deformation increases with increasing rate. Although the differences by means of maximum cortex strains are minimal, a rate dependent or stress-based failure criterion for the tau protein (or of other filaments) could be more appropriate and might alter this non-intuitive relation between poration and strain rate. Alternatively, a surface tension-based poration criterion could be investigated.
Local deformations of the axonal cortex were used as input to MD simulations of the lipid bilayer to establish model-informed axonal injury thresholds. Combining the axonal and membrane level we show that poration occurs in a range of axonal strains going from 10 to 15%. These values are in line with previously established axonal injury thresholds (18, 68–72, 105, 106). Closest agreement is found between our thresholds and those obtained with uniaxially oriented neuronal primary neuronal cultures and a spinal nerve root stretch injury mode which also observed injurious changes above 10% applied strain (105, 106). Our results are also in accordance with experimental observations that observed mechanoporation only when strains higher than 10% were uniaxially applied to a non-oriented neuronal culture (107). Thresholds obtained from experimental tissue-level injury models cover a slightly higher range of axonal strains (70, 71). However, this is probably due to the influence of axonal tortuosity and the alternation between an affine/non-affine regime (99, 100). Namely, when deforming the tissue experimentally, at least initially, part of the tissue deformation goes into straightening of the axons. Only later the axons themselves sustain a deformation.
Several studies have so far tried to assess axonal injury by means of permeability (6, 9, 10, 12, 16, 108). Cell-impermeant macromolecules, such as horseradish peroxidase or Lucifer Yellow, are commonly injected in the extracellular space, prior to the application of a mechanical insult to cell culture or the organ itself. The presence of these molecules in the intra-cellular space post-injury has been proposed as an indicator of pore formation in the axolemma. These changes are studied as a function of the magnitude of the applied strain and strain rate. However, most often results are representative of the entire culture (consisting of non-aligned axons) and not of the single cell. Hence, they cannot be compared with our results that instead relate axonal strains to permeability changes in a localized membrane volume. In a study by Kilinc et al. (10), however, permeability was observed in individual axons. In particular, membrane permeability in injured axons was found to be twice as high as in control axons. This is in agreement with our results showing that changes in water permeability can be expected to be higher than 40% (51% in presence of embedded proteins) at the onset of poration. Nevertheless, our results show that changes in permeability can be expected even below the poration threshold, namely at axonal strains lower than 10–15%. Interestingly, in a previous study by Yuen et al. (72), which utilized an aligned cell culture system to investigate mild TBI, pathological axonal alterations were observed at strains starting from εaxon =5%—a value that is more than twice as low as established axonal injury thresholds (18, 68–71). Based on our molecular simulations insight, this might be indeed justified by alterations in axolemma permeability pre-poration.
At this stage, we would like to stress that it is only the combination between these two models (the axonal and molecular-based membrane model) that allowed us to achieve such a good agreement with experimental injury thresholds. In particular, ignoring the composite nature of the axonal model, one would be tempted to assume that the strains sustained by the membrane are equivalent to the tensile strains sustained by the axon. In this scenario, according to our results (Figure 6), poration would occur at εx = εaxon ϵ [0.34, 0.51]. These strain values are considerably higher than those so far proposed for axonal injury. Our modeling effort was first aimed at an accurate representation of the axon passive mechanical response to provide axon-specific boundary conditions to simulate molecular-level membrane simulation. Secondly, an accurate representation of the membrane lipid content was deemed necessary to assess features such as membrane permeability and poration given that lipid content is known to have an influence on membrane's mechanical properties (28, 30, 31, 33). Notably, considering the strain concentrations revealed by the axonal model with a different lipid bilayer model would have yielded different results. Pure phospholipidic bilayers, in fact, have been reported not to rupture until 0.68 von Mises strain (29).
Limitations
Although our modeling approach brings new insights into the axonal injury mechanism, here we want to address some limitations. Firstly, it must be noted that no active molecular mechanism was considered in our axonal injury simulations. This was deemed a fair approximation given that the aim of this work is to assess injury events that happen in fractions of seconds, whereas biomolecular mechanisms, such as polarization/depolarization, neurofilaments transport, etc. live in the seconds-time scale. In light of this, our results are to be seen as insights into the possible mechanistic events related to primary injury, rather than damage evolution toward secondary axotomy.
It is also important to stress that our model represents a generic portion within the distal axon of an unmyelinated neuron. Several studies have so far proposed specific axonal sites such as the nodal, paranodal, internodal segments as primary sites of injury (109–111). Few studies have also highlighted the susceptibility to injury of the axon initial segment (AIS), the parasomatic region where action potentials are initiated (112–114). Due to the differences in the cytoskeleton of the AIS and the distal axon, however, our results cannot be generalized to this segment.
Microtubules failure was discarded as injury trigger in our previous publication (26). Nevertheless, primary effects on other filaments might still subsist. For example, neurofilament (NF) compaction was reported to occur as a direct effect of the mechanical insult by Meythaler et al. (115). However, a larger body of evidence suggests NFs compaction to be associated it with axolemmal permeability changes (116–118). Our model includes a representation of the NFs network, whose morphological and material characteristics are reported in the Supplementary Material. Nevertheless, we are persuaded that the current representation of this network in our FE model, despite being reasonable as a bulk, cannot reveal local phenomena related to injury. The lack of experimental information in fact prevented us from assigning specific properties to filament backbone and side-arms as well as failure properties for their connections. Hence, a priori, we cannot exclude that damage to this component of the model takes place.
At this point, it is to note that our model accounts for the effect of an embedded protein on the membrane structural features, but does not account for possible protein structural change at different strains, since the protein is described as a semi-rigid body. Nav1.1 was chosen as a representative example of axolemmal sodium channels since only this channel type's three-dimensional structure was available when this study was conducted. Nav1.1 is however, characteristic of nodal portions of myelinates axons, while Nav1.2 is found in the membrane of unmyelinated axons (78). Nevertheless, being the protein described using a semi-rigid coarse-grained model, and being Nav1.1 and Nav1.2 similar in their general structural features, the use of Nav1.1 in place of Nav1.2 should not dramatically affect the results at this level of approximation. Moreover, our protein-membrane system does not account for the inhomogeneity of proteins' distribution along the axolemma. Future studies could aim at assessing the potential effect of protein deformability, compositions or distributions on axolemma stretch susceptibility.
Moreover, the molecular systems used in this study do not account for the recruitment of lipids, which has previously been discussed as a strategy put in place by the neurons to reduce the build-up in membrane tension (34, 119). Such a mechanism was however observed in experiments applying a “gentle” osmotic perturbation rather than very fast deformations as those simulated in this study. Hence, although such a mechanism should be taken into account when studying the evolution of axonal injury, it is here deemed not to affect our considerations.
Further method specific limitations have been previously addressed for both the axonal model and for molecular simulations of membranes (26, 120). What we would like to address at this stage are possible alternatives to the methodology that was presented here to bridge the axon finite element and membrane molecular descriptions. An alternative to our approach could have been to describe also the lipid bilayer at the continuum (axonal) level and derive tensions (pressures) to be subsequently input into the molecular system. However, extracting continuum properties from molecular systems is a non-trivial problem and underlies several assumptions (64, 121–124). Equating the strains at the two different length scales, on the contrary, allowed us to prescind from a continuum description of the lipid bilayer that would have been otherwise affected by several other assumptions. Another option is to directly apply strain rates (1, 10, and 40 1/s) to the molecular-based membrane model. However, at the moment, it is computationally unfeasible to obtain results at such time scale in a reasonable time frame. Applying local strain directly to the simulation allowed us to overcome the time scale issue.
The current models were used in cascade to investigate the response to a purely uniaxial stretch injury. Although this can be comparable to a controlled in vitro loading mode, it does not directly correlate with the in vivo scenario. Strains in the brain are in fact represented by a 3D tensor. Current finite element head models resolve, with slightly different methods, this tensor in the fiber direction and use this one-dimensional quantity as injury predictor (125–128). In future studies, the uniaxial condition could be dropped to study how an axonal model embedded in a matrix responds to a complete strain tensor. As a result of the different boundary conditions the axonal model and its axolemma might undergo deformations which are different from those observed in the present study.
Conclusions
In this study, we successfully applied a top-to-bottom approach to better characterize the onset of axonal injury. By bridging continuum and molecular descriptions, we have provided quantitative evidences showing that mechanoporation of the axolemma is an event that cannot be excluded in a typical axonal injury scenario. In addition, we showed that pre-poration strain levels were characterized by an increased water permeability and interdigitated states. All in all, our study results provide increased knowledge of the axonal injury mechanism and have therefore potential positive implications for the development or optimization of neuroprotective measures targeting membrane integrity in the treatment of axonal injuries and brain injuries in general.
Data Availability Statement
The datasets generated for this study are available on request to the corresponding author.
Author Contributions
AM, MS, AV, and SK designed the study. AM performed and analyzed the finite element simulations, while MS and AV performed and analyzed the molecular dynamics simulation. AM, MS, and AV prepared the manuscript. All authors critically reviewed the manuscript.
Funding
The authors thank the Swedish Research Council (VR-2016-05314), the European Union Horizon 2020 Research and Innovation Framework Programme under the Marie Skłodowska-Curie (grant agreement no. 642662), and the Swedish National Infrastructure for Computing (SNIC2017-1-491 and SNIC2018-3-548) for the support.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fneur.2020.00025/full#supplementary-material
References
1. Menon DK, Schwab K, Wright DW, Maas AI. Position statement: definition of traumatic brain injury. Arch Phys Med Rehabil. (2010) 91:1637–40. doi: 10.1016/j.apmr.2010.05.017
2. Taylor CA, Bell JM, Breiding MJ, Xu L. Traumatic brain injury-related emergency department visits, hospitalizations, and deaths - United States, 2007 and 2013. MMWR Surveill Summ. (2017) 66:1–16. doi: 10.15585/mmwr.ss6609a1
3. Peeters W, van den Brande R, Polinder S, Brazinova A, Steyerberg EW, Lingsma HF, et al. Epidemiology of traumatic brain injury in Europe. Acta Neurochir. (2015) 157:1683–96. doi: 10.1007/s00701-015-2512-7
4. Johnson VE, Stewart W, Smith DH. Axonal pathology in traumatic brain injury. Exp Neurol. (2013) 246:35–43. doi: 10.1016/j.expneurol.2012.01.013
5. Hill CS, Coleman MP, Menon DK. Traumatic axonal injury: mechanisms and translational opportunities. Trends Neurosci. (2016) 39:311–24. doi: 10.1016/j.tins.2016.03.002
6. Pettus EH, Christman CW, Giebel ML, Povlishock JT. Traumatically induced altered membrane permeability: its relationship to traumatically induced reactive axonal change. J Neurotrauma. (1994) 11:507–22. doi: 10.1089/neu.1994.11.507
7. LaPlaca MC, Thibault LE. An in vitro traumatic injury model to examine the response of neurons to a hydrodynamically-induced deformation. Ann Biomed Eng. (1997) 25:665–77. doi: 10.1007/BF02684844
8. Fitzpatrick MO, Maxwell WL, Graham DI. The role of the axolemma in the initiation of traumatically induced axonal injury. J Neurol Neurosurg Psychiatr. (1998) 64:285–7. doi: 10.1136/jnnp.64.3.285
9. Geddes DM, Cargill RS, LaPlaca MC. Mechanical stretch to neurons results in a strain rate and magnitude-dependent increase in plasma membrane permeability. J Neurotrauma. (2003) 20:1039–49. doi: 10.1089/089771503770195885
10. Kilinc D, Gallo G, Barbee KA. Mechanically-induced membrane poration causes axonal beading and localized cytoskeletal damage. Exp Neurol. (2008) 212:422–30. doi: 10.1016/j.expneurol.2008.04.025
11. Cullen DK, LaPlaca MC. Neuronal response to high rate shear deformation depends on heterogeneity of the local strain field. J Neurotrauma. (2006) 23:1304–19. doi: 10.1089/neu.2006.23.1304
12. Cullen DK, Vernekar VN, LaPlaca MC. Trauma-induced plasmalemma disruptions in three-dimensional neural cultures are dependent on strain modality and rate. J Neurotrauma. (2011) 28:2219–33. doi: 10.1089/neu.2011.1841
13. LaPlaca MC, Lessing MC, Prado GR, Zhou R, Tate CC, Geddes-Klein D, et al. Mechanoporation is a potential indicator of tissue strain and subsequent degeneration following experimental traumatic brain injury. Clin Biomech. (2019) 64:2–13. doi: 10.1016/j.clinbiomech.2018.05.016
14. LaPlaca MC, Prado GR. Neural mechanobiology and neuronal vulnerability to traumatic loading. J Biomech. (2010) 43:71–8. doi: 10.1016/j.jbiomech.2009.09.011
15. Wolf JA, Stys PK, Lusardi T, Meaney D, Smith DH. Traumatic axonal injury induces calcium influx modulated by tetrodotoxin-sensitive sodium channels. J Neurosci. (2001) 21:1923–30. doi: 10.1523/JNEUROSCI.21-06-01923.2001
16. Stone JR, Okonkwo DO, Dialo AO, Rubin DG, Mutlu LK, Povlishock JT, et al. Impaired axonal transport and altered axolemmal permeability occur in distinct populations of damaged axons following traumatic brain injury. Exp Neurol. (2004) 190:59–69. doi: 10.1016/j.expneurol.2004.05.022
17. Farkas O. Mechanoporation Induced by diffuse traumatic brain injury: an irreversible or reversible response to injury? J Neurosci. (2006) 26:3130–40. doi: 10.1523/JNEUROSCI.5119-05.2006
18. Tang-Schomer MD, Patel AR, Baas PW, Smith DH. Mechanical breaking of microtubules in axons during dynamic stretch injury underlies delayed elasticity, microtubule disassembly, and axon degeneration. FASEB J. (2010) 24:1401–10. doi: 10.1096/fj.09-142844
19. Tang-Schomer MD, Johnson VE, Baas PW, Stewart W, Smith DH. Partial interruption of axonal transport due to microtubule breakage accounts for the formation of periodic varicosities after traumatic axonal injury. Exp Neurol. (2012) 233:364–72. doi: 10.1016/j.expneurol.2011.10.030
20. Peter SJ, Mofrad MRK. Computational modeling of axonal microtubule bundles under tension. Biophys J. (2012) 102:749–57. doi: 10.1016/j.bpj.2011.11.4024
21. Soheilypour M, Peyro M, Peter SJ, Mofrad MRK. Buckling behavior of individual and bundled microtubules. Biophys J. (2015) 108:1718–26. doi: 10.1016/j.bpj.2015.01.030
22. Lazarus C, Soheilypour M, Mofrad MRK. Torsional behavior of axonal microtubule bundles. Biophys J. (2015) 109:231–9. doi: 10.1016/j.bpj.2015.06.029
23. Ahmadzadeh H, Smith DH, Shenoy VB. Viscoelasticity of tau proteins leads to strain rate-dependent breaking of microtubules during axonal stretch injury: predictions from a mathematical model. Biophys J. (2014) 106:1123–33. doi: 10.1016/j.bpj.2014.01.024
24. Ahmadzadeh H, Smith DH, Shenoy VB. Mechanical effects of dynamic binding between tau proteins on microtubules during axonal injury. Biophys J. (2015) 109:2328–37. doi: 10.1016/j.bpj.2015.09.010
25. de Rooij R, Kuhl E. Physical biology of axonal damage. Front Cell Neurosci. (2018) 12:144. doi: 10.3389/fncel.2018.00144
26. Montanino A, Kleiven S. Utilizing a structural mechanics approach to assess the primary effects of injury loads onto the axon and its components. Front Neurol. (2018) 9:643. doi: 10.3389/fneur.2018.00643
27. Mohandas N, Evans E. Mechanical properties of the red cell membrane in relation to molecular structure and genetic defects. Annu Rev Biophys Biomol Struct. (1994) 23:787−818. doi: 10.1146/annurev.bb.23.060194.004035
28. Shigematsu T, Koshiyama K, Wada S. Effects of stretching speed on mechanical rupture of phospholipid/cholesterol bilayers: molecular dynamics simulation. Sci Rep. (2015) 5:15369. doi: 10.1038/srep15369
29. Murphy MA, Horstemeyer MF, Gwaltney SR, Stone T, LaPlaca M, Liao J, et al. Nanomechanics of phospholipid bilayer failure under strip biaxial stretching using molecular dynamics. Model Simul Mater Sci Eng. (2016) 24:055008. doi: 10.1088/0965-0393/24/5/055008
30. Saeedimasine M, Montanino A, Kleiven S, Villa A. Role of lipid composition on the structural and mechanical features of axonal membranes: a molecular simulation study. Sci Rep. (2019) 9:8000. doi: 10.1038/s41598-019-44318-9
31. Harayama T, Riezman H. Understanding the diversity of membrane lipid composition. Nat Rev Mol Cell Biol. (2018) 19:281. doi: 10.1038/nrm.2017.138
32. Evans EA, Waugh R, Melnik L. Elastic area compressibility modulus of red cell membrane. Biophys J. (1976) 16:585–95. doi: 10.1016/S0006-3495(76)85713-X
33. Needham D, Nunn RS. Elastic deformation and failure of lipid bilayer membranes containing cholesterol. Biophys J. (1990) 58:997–1009. doi: 10.1016/S0006-3495(90)82444-9
34. Pullarkat PA, Dommersnes P, Fernández P, Joanny J-F, Ott A. Osmotically driven shape transformations in axons. Phys Rev Lett. (2006) 96:48104. doi: 10.1103/PhysRevLett.96.048104
35. Wann KT. Neuronal sodium and potassium channels: structure and function. Br J Anaesth. (1993) 71:2–14. doi: 10.1093/bja/71.1.2
36. Hu H, Jonas P. A supercritical density of Na+ channels ensures fast signaling in GABAergic interneuron axons. Nat Neurosci. (2014) 17:686. doi: 10.1038/nn.3678
37. Iwata A. Traumatic axonal injury induces proteolytic cleavage of the voltage-gated sodium channels modulated by tetrodotoxin and protease inhibitors. J Neurosci. (2004) 24:4605–13. doi: 10.1523/JNEUROSCI.0515-03.2004
38. Zhang L, Zhang Z, Jasa J, Li D, Cleveland RO, Negahban M, et al. Molecular dynamics simulations of heterogeneous cell membranes in response to uniaxial membrane stretches at high loading rates. Sci Rep. (2017) 7:8316. doi: 10.1038/s41598-017-06827-3
39. Reeves TM, Phillips LL, Povlishock JT. Myelinated and unmyelinated axons of the corpus callosum differ in vulnerability and functional recovery following traumatic brain injury. Exp Neurol. (2005) 196:126–37. doi: 10.1016/j.expneurol.2005.07.014
40. Staal JA, Vickers JC. Selective vulnerability of non-myelinated axons to stretch injury in an in vitro co-culture system. J Neurotrauma. (2011) 28:841–7. doi: 10.1089/neu.2010.1658
41. Karcher H, Lammerding J, Huang H, Lee RT, Kamm RD, Kaazempur-Mofrad MR. A three-dimensional viscoelastic model for cell deformation with experimental verification. Biophys J. (2003) 85:3336–49. doi: 10.1016/S0006-3495(03)74753-5
42. Dao M, Lim CT, Suresh S. Mechanics of the human red blood cell deformed by optical tweezers. J Mech Phys Solids. (2003) 51:2259–80. doi: 10.1016/j.jmps.2003.09.019
43. Barreto S, Clausen CH, Perrault CM, Fletcher DA, Lacroix D. A multi-structural single cell model of force-induced interactions of cytoskeletal components. Biomaterials. (2013) 34:6119–26. doi: 10.1016/j.biomaterials.2013.04.022
44. Bavi N, Nakayama Y, Bavi O, Cox CD, Qin Q-H, Martinac B. Biophysical implications of lipid bilayer rheometry for mechanosensitive channels. Proc Natl Acad Sci USA. (2014) 111:13864–9. doi: 10.1073/pnas.1409011111
45. Jérusalem A, Dao M. Continuum modeling of a neuronal cell under blast loading. Acta Biomater. (2012) 8:3360–71. doi: 10.1016/j.actbio.2012.04.039
46. Bathe M, Heussinger C, Claessens MMAE, Bausch AR, Frey E. Cytoskeletal bundle mechanics. Biophys J. (2008) 94:2955–64. doi: 10.1529/biophysj.107.119743
47. Chen X, Cui Q, Tang Y, Yoo J, Yethiraj A. Gating mechanisms of mechanosensitive channels of large conductance, i: a continuum mechanics-based hierarchical framework. Biophys J. (2008) 95:563–80. doi: 10.1529/biophysj.107.128488
48. Zhu F, Gatti DL, Yang KH. Nodal versus total axonal strain and the role of cholesterol in traumatic brain injury. J Neurotrauma. (2016) 33:859–70. doi: 10.1089/neu.2015.4007
49. Cinelli I, Destrade M, Duffy M, McHugh P. Electrothermal equivalent three-dimensional finite-element model of a single neuron. IEEE Trans Biomed Eng. (2018) 65:1373–81. doi: 10.1109/TBME.2017.2752258
50. Kwong MT, Bianchi F, Malboubi M, García-Grajales JA, Homsi L, Thompson M, et al. 3D finite element formulation for mechanical–electrophysiological coupling in axonopathy. Comput Methods Appl Mech Eng. (2019) 346:1025–50. doi: 10.1016/j.cma.2018.09.006
51. Bernal R, Pullarkat PA, Melo F. Mechanical properties of axons. Phys Rev Lett. (2007) 99:18301. doi: 10.1103/PhysRevLett.99.018301
52. Dennerll TJ, Joshi HC, Steel VL, Buxbaum RE, Heidemann SR. Tension and compression in the cytoskeleton of PC-12 neurites II: quantitative measurements. J Cell Biol. (1988) 107:665–74. doi: 10.1083/jcb.107.2.665
53. Rajagopalan J, Tofangchi A, Saif MTA. Drosophila neurons actively regulate axonal tension in vivo. Biophys J. (2010) 99:3208–15. doi: 10.1016/j.bpj.2010.09.029
54. Yu W, Baas PW. Changes in microtubule number and length during axon differentiation. J Neurosci. (1994) 14:2818–29. doi: 10.1523/JNEUROSCI.14-05-02818.1994
55. Bray D, Bunge MB. Serial analysis of microtubules in cultured rat sensory axons. J Neurocytol. (1981) 10:589–605. doi: 10.1007/BF01262592
56. Fadić R, Vergara J, Alvarez J. Microtubules and caliber of central and peripheral processes of sensory axons. J. Comp. Neurol. (1985) 236:258–64. doi: 10.1002/cne.902360209
57. Malbouisson AMB, Ghabriel MN, Allt G. Axonal microtubules: a computer-linked quantitative analysis. Anat Embryol. (1985) 171:339–44. doi: 10.1007/BF00347022
58. Leterrier C, Dubey P, Roy S. The nano-architecture of the axonal cytoskeleton. Nat Rev Neurosci. (2017) 18:713. doi: 10.1038/nrn.2017.129
59. Evans EA, Hochmuth RM. A solid-liquid composite model of the red cell membrane. J Membr Biol. (1976) 30:351–62. doi: 10.1007/BF01869676
60. Xu K, Zhong G, Zhuang X. Actin, spectrin, and associated proteins form a periodic cytoskeletal structure in axons. Science. (2013) 339:452–6. doi: 10.1126/science.1232251
61. Datar A, Bornschlögl T, Bassereau P, Prost J, Pullarkat PA. Dynamics of membrane tethers reveal novel aspects of cytoskeleton-membrane interactions in axons. Biophys J. (2015) 108:489–97. doi: 10.1016/j.bpj.2014.11.3480
62. Zhang Y, Abiraman K, Li H, Pierce DM, Tzingounis AV, Lykotrafitis G. Modeling of the axon membrane skeleton structure and implications for its mechanical properties. PLoS Comput Biol. (2017) 13:1–22. doi: 10.1371/journal.pcbi.1005407
63. Pullarkat PA, Fernández PA, Ott A. Rheological properties of the Eukaryotic cell cytoskeleton. Phys Rep. (2007) 449:29–53. doi: 10.1016/j.physrep.2007.03.002
64. Zhang L, Jasa J, Gazonas G, Jérusalem A, Negahban M. Extracting continuum-like deformation and stress from molecular dynamics simulations. Comput Methods Appl Mech Eng. (2015) 283:1010–31. doi: 10.1016/j.cma.2014.10.018
65. Duke T, Leibler S. Motor protein mechanics: a stochastic model with minimal mechanochemical coupling. Biophys J. (1996) 71:1235–47. doi: 10.1016/S0006-3495(96)79323-2
66. Bain AC, Raghupathi R, Meaney DF. Dynamic stretch correlates to both morphological abnormalities and electrophysiological impairment in a model of traumatic axonal injury. J Neurotrauma. (2001) 18:499–511. doi: 10.1089/089771501300227305
67. Geddes-Klein DM, Schiffman KB, Meaney DF. Mechanisms and consequences of neuronal stretch injury in vitro differ with the model of trauma. J Neurotrauma. (2006) 23:193–204. doi: 10.1089/neu.2006.23.193
68. Margulies SS, Thibault LE. A proposed tolerance criterion for diffuse axonal injury in man. J Biomech. (1992) 25:917–23. doi: 10.1016/0021-9290(92)90231-O
69. Gennarelli TA, Thibault LE, Graham DI. Diffuse axonal injury: an important form of traumatic brain damage. Neurosci. (1998) 4:202–15. doi: 10.1177/107385849800400316
70. Bain AC, Meaney DF. Tissue-level thresholds for axonal damage in an experimental model of central nervous system white matter injury. J Biomech Eng. (2000) 122:615–22. doi: 10.1115/1.1324667
71. Elkin BS, Morrison B. Region-specific tolerance criteria for the living brain. In: 51st Stapp Car Crash Conference. San Diego, CA: The Stapp Association (2007). doi: 10.4271/2007-22-0005
72. Yuen TJ, Browne KD, Iwata A, Smith DH. Sodium channelopathy induced by mild axonal trauma worsens outcome after a repeat injury. J Neurosci Res. (2009) 87:3620–5. doi: 10.1002/jnr.22161
73. Marrink SJ, de Vries AH, Mark AE. Coarse grained model for semiquantitative lipid simulations. J Phys Chem B. (2004) 108:750–60. doi: 10.1021/jp036508g
74. Marrink SJ, Risselada HJ, Yefimov S, Tieleman DP, De Vries AH. The MARTINI force field: coarse grained model for biomolecular simulations. J Phys Chem B. (2007) 111:7812–24. doi: 10.1021/jp071097f
75. Yesylevskyy SO, Schäfer LV, Sengupta D, Marrink SJ. Polarizable water model for the coarse-grained MARTINI force field. PLoS Comput Biol. (2010) 6:e1000810. doi: 10.1371/journal.pcbi.1000810
76. Ingólfsson HI, Melo MN, Van Eerden FJ, Arnarez C, Lopez CA, Wassenaar TA, et al. Lipid organization of the plasma membrane. J Am Chem Soc. (2014) 136:14554–9. doi: 10.1021/ja507832e
77. Catterall WA, Goldin AL, Waxman SG. International union of pharmacology. XLVII nomenclature and structure-function relationships of voltage-gated sodium channels. Pharmacol Rev. (2005) 57:397–409. doi: 10.1124/pr.57.4.4
78. Westenbroek RE, Merrick DK, Catterall WA. Differential subcellular localization of the RI and RII Na+ channel subtypes in central neurons. Neuron. (1989) 3:695–704. doi: 10.1016/0896-6273(89)90238-9
79. Duflocq A, Le Bras B, Bullier E, Couraud F, Davenne M. Nav1. 1 is predominantly expressed in nodes of Ranvier and axon initial segments. Mol Cell Neurosci. (2008) 39:180–92. doi: 10.1016/j.mcn.2008.06.008
80. Malo MS, Blanchard BJ, Andresen JM, Srivastava K, Chen X-N, Li X, et al. Localization of a putative human brain sodium channel gene (SCN1A) to chromosome band 2q24. Cytogenet Genome Res. (1994) 67:178–86. doi: 10.1159/000133818
81. Kelley LA, Mezulis S, Yates CM, Wass MN, Sternberg MJE. The Phyre2 web portal for protein modeling, prediction and analysis. Nat Protoc. (2015) 10:845. doi: 10.1038/nprot.2015.053
82. Shen H, Zhou Q, Pan X, Li Z, Wu J, Yan N. Structure of a eukaryotic voltage-gated sodium channel at near-atomic resolution. Science. (2017) 355:eaal4326. doi: 10.1126/science.aal4326
83. Eswar N, Webb B, Marti-Renom MA, Madhusudhan MS, Eramian D, Shen M, et al. Comparative protein structure modeling using modeller. Curr Protoc Bioinforma. (2006) 15:5.6.1–30. doi: 10.1002/0471250953.bi0506s15
84. Lomize MA, Pogozheva ID, Joo H, Mosberg HI, Lomize AL. OPM database and PPM web server: resources for positioning of proteins in membranes. Nucleic Acids Res. (2011) 40:D370–76. doi: 10.1093/nar/gkr703
85. Pastor RW, Mackerell AD. Development of the CHARMM Force Field for Lipids. J Phys Chem Lett. (2011) 2:1526–32. doi: 10.1021/jz200167q
86. Monticelli L, Kandasamy SK, Periole X, Larson RG, Tieleman DP, Marrink S-J. The MARTINI coarse-grained force field: extension to proteins. J Chem Theory Comput. (2008) 4:819–34. doi: 10.1021/ct700324x
87. Periole X, Cavalli M, Marrink S-J, Ceruso MA. Combining an elastic network with a coarse-grained molecular force field: structure, dynamics, and intermolecular recognition. J Chem Theory Comput. (2009) 5:2531–43. doi: 10.1021/ct9002114
88. Marrink SJ, de Vries AH, Tieleman DP. Lipids on the move: simulations of membrane pores, domains, stalks and curves. Biochim Biophys Acta - Biomembr. (2009) 1788:149–68. doi: 10.1016/j.bbamem.2008.10.006
89. Abraham MJ, Murtola T, Schulz R, Páll S, Smith JC, Hess B, et al. GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX. (2015) 1:19–25. doi: 10.1016/j.softx.2015.06.001
90. Bussi G, Zykova-Timan T, Parrinello M. Isothermal-isobaric molecular dynamics using stochastic velocity rescaling. J Chem Phys. (2009) 130:74101. doi: 10.1063/1.3073889
91. Parrinello M, Rahman A. Polymorphic transitions in single crystals: a new molecular dynamics method. J Appl Phys. (1981) 52:7182–90. doi: 10.1063/1.328693
92. Páll S, Hess B. A flexible algorithm for calculating pair interactions on SIMD architectures. Comput Phys Commun. (2013) 184:2641–50. doi: 10.1016/j.cpc.2013.06.003
93. Tironi IG, Sperb R, Smith PE, van Gunsteren WF. A generalized reaction field method for molecular dynamics simulations. J Chem Phys. (1995) 102:5451–9. doi: 10.1063/1.469273
94. Humphrey W, Dalke A, Schulten K. VMD: visual molecular dynamics. J Mol Graph. (1996) 14:33–8. doi: 10.1016/0263-7855(96)00018-5
95. Shinoda W. Permeability across lipid membranes. Biochim Biophys Acta - Biomembr. (2016) 1858:2254–65. doi: 10.1016/j.bbamem.2016.03.032
96. Walter A, Gutknecht J. Permeability of small nonelectrolytes through lipid bilayer membranes. J Membr Biol. (1986) 90:207–17. doi: 10.1007/BF01870127
97. Overton E. Ueber die osmotischen eigenschaften der zelle in ihrer bedeutung für die toxikologie und pharmakologie. Zeitschrift für Phys Chemie. (1897) 22:189–209. doi: 10.1515/zpch-1897-2220
98. Waugh R, Evans EA. Thermoelasticity of red blood cell membrane. Biophys J. (1979) 26:115–31. doi: 10.1016/S0006-3495(79)85239-X
99. Bain AC, Shreiber DI, Meaney DF. Modeling of microstructural kinematics during simple elongation of central nervous system tissue. J Biomech Eng. (2004) 125:798–804. doi: 10.1115/1.1632627
100. Singh S, Pelegri AA, Shreiber DI. Characterization of the three-dimensional kinematic behavior of axons in central nervous system white matter. Biomech Model Mechanobiol. (2015) 14:1303–15. doi: 10.1007/s10237-015-0675-z
101. Pan Y, Sullivan D, Shreiber D, Pelegri A. Finite element modeling of CNS white matter kinematics: use of a 3D RVE to determine material properties. Front Bioeng Biotechnol. (2013) 1:19. doi: 10.3389/fbioe.2013.00019
102. Chetta J, Kye C, Shah SB. Cytoskeletal dynamics in response to tensile loading of mammalian axons. Cytoskeleton. (2010) 67:650–65. doi: 10.1002/cm.20478
103. de Rooij R, Kuhl E. Microtubule polymerization and cross-link dynamics explain axonal stiffness and damage. Biophys J. (2018) 114:201–12. doi: 10.1016/j.bpj.2017.11.010
104. Bianchi F, Pereno V, George JH, Thompson MS, Ye H. Membrane mechanical properties regulate the effect of strain on spontaneous electrophysiology in human iPSC-derived. Neurons. (2019) 404:165–74. doi: 10.1016/j.neuroscience.2019.02.014
105. Singh A, Kallakuri S, Chen C, Cavanaugh JM. Structural and functional changes in nerve roots due to tension at various strains and strain rates: an in-vivo study. J Neurotrauma. (2009) 26:627–40. doi: 10.1089/neu.2008.0621
106. Nakadate H, Kurtoglu E, Furukawa H, Oikawa S, Aomura S, Kakuta A, et al. Strain-rate Dependency of axonal tolerance for uniaxial stretching. In: 61st Stapp Car Crash Conference. Charleston, SC: The Stapp Association. (2017). doi: 10.4271/2017-22-0003
107. Hemphill MA, Dabiri BE, Gabriele S, Kerscher L, Franck C, Goss JA, et al. A possible role for integrin signaling in diffuse axonal injury. PLoS ONE. (2011) 6:e22899. doi: 10.1371/journal.pone.0022899
108. Shi R, Pryor JD. Pathological changes of isolated spinal cord axons in response to mechanical stretch. Neuroscience. (2002) 110:765–77. doi: 10.1016/S0306-4522(01)00596-6
109. Maxwell WL. Histopathological changes at central nodes of ranvier after stretch-injury. Microsc. Res. Tech. (1996) 34:522–35. doi: 10.1002/(SICI)1097-0029(19960815)34:6<522::AID-JEMT4>3.0.CO;2-L
110. Pettus EH, Povlishock JT. Characterization of a distinct set of intra-axonal ultrastructural changes associated with traumatically induced alteration in axolemmal permeability. Brain Res. (1996) 722:1–11. doi: 10.1016/0006-8993(96)00113-8
111. Reeves TM, Greer JE, Vanderveer AS, Phillips LL. Proteolysis of submembrane cytoskeletal proteins ankyrin-G and αII-spectrin following diffuse brain injury: a role in white matter vulnerability at nodes of ranvier. Brain Pathol. (2010) 20:1055–68. doi: 10.1111/j.1750-3639.2010.00412.x
112. Schafer DP, Jha S, Liu F, Akella T, McCullough LD, Rasband MN. Disruption of the axon initial segment cytoskeleton is a new mechanism for neuronal injury. J Neurosci. (2009) 29:13242–54. doi: 10.1523/JNEUROSCI.3376-09.2009
113. Buffington SA, Rasband MN. The axon initial segment in nervous system disease and injury. Eur J Neurosci. (2011) 34:1609–19. doi: 10.1111/j.1460-9568.2011.07875.x
114. Greer JE, Hånell A, McGinn MJ, Povlishock JT. Mild traumatic brain injury in the mouse induces axotomy primarily within the axon initial segment. Acta Neuropathol. (2013) 126:59–74. doi: 10.1007/s00401-013-1119-4
115. Meythaler JM, Peduzzi JD, Eleftheriou E, Novack TA. Current concepts: diffuse axonal injury–associated traumatic brain injury. Arch Phys Med Rehabil. (2001) 82:1461–71. doi: 10.1053/apmr.2001.25137
116. Povlishock JT, Marmarou A, McIntosh T, Trojanowski JQ, Moroi J. Impact acceleration injury in the rat: evidence for focal axolemmal change and related neurofilament sidearm alteration. J Neuropathol Exp Neurol. (1997) 56:347–59. doi: 10.1097/00005072-199704000-00003
117. Okonkwo DO, Pettus EH, Moroi J, Povlishock JT. Alteration of the neurofilament sidearm and its relation to neurofilament compaction occurring with traumatic axonal injury. Brain Res. (1998) 784:1–6. doi: 10.1016/S0006-8993(97)01075-5
118. Marmarou CR, Povlishock JT. Administration of the immunophilin ligand FK506 differentially attenuates neurofilament compaction and impaired axonal transport in injured axons following diffuse traumatic brain injury. Exp Neurol. (2006) 197:353–62. doi: 10.1016/j.expneurol.2005.10.003
119. Dai J, Sheetz MP, Wan X, Morris CE. Membrane tension in swelling and shrinking molluscan neurons. J Neurosci. (1998) 18:6681–92. doi: 10.1523/JNEUROSCI.18-17-06681.1998
120. Marrink SJ, Corradi V, Souza PCT, Ingólfsson HI, Tieleman DP, Sansom MSP. Computational modeling of realistic cell membranes. Chem. Rev. (2019) 119:6184–226. doi: 10.1021/acs.chemrev.8b00460
121. Subramaniyan AK, Sun CT. Continuum interpretation of virial stress in molecular simulations. Int J Solids Struct. (2008) 45:4340–6. doi: 10.1016/j.ijsolstr.2008.03.016
122. Zhou M. A new look at the atomic level virial stress: on continuum-molecular system equivalence. Proc R Soc London Ser A Math Phys Eng Sci. (2003) 459:2347–92. doi: 10.1098/rspa.2003.1127
123. Zimmerman JA, Webb EB, Hoyt JJ, Jones RE, Klein PA, Bammann DJ. - Evaluation of continuum stress in atomistic simulation. In: Bathe KJ, editor. Computational Fluid and Solid Mechanics 2003. Cambridge, MA: Elsevier Science Ltd (2003). p. 804–7. doi: 10.1016/B978-008044046-0.50196-2
124. Vanegas JM, Torres-Sánchez A, Arroyo M. Importance of force decomposition for local stress calculations in biomembrane molecular simulations. J Chem Theory Comput. (2014) 10:691–702. doi: 10.1021/ct4008926
125. Wright RM, Ramesh KT. An axonal strain injury criterion for traumatic brain injury. Biomech Model Mechanobiol. (2012) 11:245–60. doi: 10.1007/s10237-011-0307-1
126. Giordano C, Kleiven S. Evaluation of axonal strain as a predictor for mild traumatic brain injuries using finite element modeling. In: 58th Stapp Car Crash Conference. The Stapp Association. (2014). doi: 10.4271/2014-22-0002
127. Ji S, Zhao W, Ford JC, Beckwith JG, Bolander RP, Greenwald RM, et al. Group-wise evaluation and comparison of white matter fiber strain and maximum principal strain in sports-related concussion. J Neurotrauma. (2015) 32:441–54. doi: 10.1089/neu.2013.3268
Keywords: mechanoporation, axolemma, axonal injury, membrane permeability, traumatic brain injury, finite element
Citation: Montanino A, Saeedimasine M, Villa A and Kleiven S (2020) Localized Axolemma Deformations Suggest Mechanoporation as Axonal Injury Trigger. Front. Neurol. 11:25. doi: 10.3389/fneur.2020.00025
Received: 23 October 2019; Accepted: 09 January 2020;
Published: 30 January 2020.
Edited by:
Mattias K. Sköld, Uppsala University, SwedenReviewed by:
James C. Vickers, University of Tasmania, AustraliaKeiichiro Susuki, Wright State University, United States
Copyright © 2020 Montanino, Saeedimasine, Villa and Kleiven. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Annaclaudia Montanino, annacl2@kth.se