Elucidating Axonal Injuries Through Molecular Modelling of Myelin Sheaths and Nodes of Ranvier.

Around half of the traumatic brain injuries are thought to be axonal damage. Disruption of the cellular membranes, or alternatively cytoskeletal damage has been suggested as possible injury trigger. Here, we have used molecular models to have a better insight on the structural and mechanical properties of axon sub-cellular components. We modelled myelin sheath and node of Ranvier as lipid bilayers at a coarse grained level. We built ex-novo a model for the myelin. Lipid composition and lipid saturation were based on the available experimental data. The model contains 17 different types of lipids, distributed asymmetrically between two leaflets. Molecular dynamics simulations were performed to characterize the myelin and node-of-Ranvier bilayers at equilibrium and under deformation and compared to previous axolemma simulations. We found that the myelin bilayer has a slightly higher area compressibility modulus and higher rupture strain than node of Ranvier. Compared to the axolemma in unmyelinated axon, mechanoporation occurs at 50% higher strain in the myelin and at 23% lower strain in the node of Ranvier in myelinated axon. Combining the results with finite element simulations of the axon, we hypothesizes that myelin does not rupture at the thresholds proposed in the literature for axonal injury while rupture may occur at the node of Ranvier. The findings contribute to increases our knowledge of axonal sub-cellular components and help to understand better the mechanism behind axonal brain injury.

1 Introduction

Traumatic brain injury (TBI) (Menon et al., 2010) is an injury to the brain caused by an external force. Causes include falls, vehicle collisions and violence. One of the most common consequence of a TBI is diffuse axonal injury, a multifocal damage to white matter (Johnson et al., 2013). Such an injury is invisible to conventional brain imaging, and can only be histologically diagnosed and one of its hallmarks is the presence of axonal swellings (Hill et al., 2016).

Axons are long projection of the nerve cell surrounded by a membrane, called axolemma. Axolemma separates the interior of the axon from the outside environments. In the nervous system, axons may be myelinated, or unmyelinated. When the axon is myelinated, a extra multilamellar membrane, called myelin sheath, insulates segments of axon. The sheath consists of repeating units of double bilayers separated by 3–4 nm-thick aqueous layers that alternate between the cytoplasmic and extracellular faces of membranes (Inouye and Kirschner, 1988a). Dehydrated myelin is unusual in that it is composed of 75–80% lipid and 20–25% protein by weight, compared with around 50% of lipids in most other cell membranes (Williams et al., 1993). Multiple lipids make up the myelin sheath, and each sheath contributes to the structure, adhesive stability, and possibly the pathogenesis of the myelin membrane. The asymmetric distribution of lipid composition on the cytoplasmic and extracellular faces likely also plays an important role (Inouye and Kirschner, 1988b). Axolemma and myelin compositions differ both in lipid type and degree of saturation. Table 1 summarizes the available experimental data on lipid composition for axolemma and myelin. The most striking feature of myelin lipid composition is the enrichment in glycolipids together with long-chain fatty acids (DeVries et al., 1981).

TABLE 1

Experimental and modelled lipid composition in weight fraction for different types of membrane.

Lipid type Membrane	PC	PE	SM	PS	Glycolipid	CHOL	Other
Myelin (exp) a	0.11	0.17	0.06	0.06	0.28	0.28	0.04
	(0.01)	(0.01)	(0.03)	(0.01)	(0.03)	(0.01)	(0.02)
Myelin model	0.12	0.22	0.03	0.09	0.23	0.29	0.02
Axolemma (exp) b	0.20	0.19	0.08	0.06	0.20	0.24	0.10
	(0.09)	(0.07)	(0.05)	(0.03)	(0.05)	(0.03)	(0.07)
Plasma model c	0.30	0.17	0.17	0.06	0.05	0.17	0.08

For composition observed experimentally the average values is reported together with the maximum deviations from the mean values (in parentheses).

data for the central nervous system of human, bovine and rat from Rasband and Macklin (2012) and references inside.

data from central and peripheral nervous system from Camejo et al. (1969); Zambrano et al. (1971); DeVries et al. (1981); DeVries et al. (1999).

plasma membrane model of Ingólfsson et al. (2014) was used to describe the axolemma.

Short unmyelinated segments, nodes of Ranvier, occur periodically between segments of the myelin sheath in myelineted axons. In nodes of Ranvier, the axolemma is directly exposed to the extracellular space. The nodes are uninsulated and highly enriched in ion channels (Westenbroek et al., 1989; Catterall et al., 2005), allowing them to participate in the exchange of ions required to regenerate the action potential. One of the voltage-gated ion channels, embedded in myelinated axons, is sodium channel protein type subunit alpha, Nav1.1 (Duflocq et al., 2008). In axonal membrane, the density of sodium channels varies between 5 and 3,000 channels/μm2: lower densities are observed in unmyelinated axons, while higher densities are found in the nodal portions of myelinated axons (Wann, 1993; Hu and Jonas, 2014). Interestingly, sodium channels have also been proposed to influence the injury response (Iwata, 2004).

Experimental observations have so far led to the formulation of two main theories regarding the cellular primary injury mechanism. Disruption of the axolemma (Pettus and Povlishock, 1996; LaPlaca and Thibault, 1997; Fitzpatrick et al., 1998), or alternatively cytoskeletal damage (Tang-Schomer et al., 2012) has been suggested mainly as injury trigger. However, using a purely mechanical approach we discarded microtubule damage as injury trigger and revealed instead high level of strains on the axonal membrane (Montanino and Kleiven, 2018). To further investigate the molecular level effects of such strains, we bridged the finite element model of the axon with a molecular-based membrane model (Montanino et al., 2020). Despite the approximation of the models, we showed 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 at which poration occurs (10–12%) agrees well with the thresholds obtained both with in-vitro (Hemphill et al., 2011; Nakadate et al., 2017) and in-vivo/ex-vivo (Bain and Meaney, 2000; Singh et al., 2015) stretch injury experiments and allows us to provide quantitative evidences that do not exclude pore formation in the membrane as a result of trauma.

When investigating the progression of axonal injury in the white matter, several studies have observed impairments at the myelin level in the form of delamination of the myelin lamellae (Maxwell, 2013) and general loss of myelin, or demyelination, (Mierzwa et al., 2015; Mu et al., 2019). While white matter degeneration is undoubtedly a distinctive feature of traumatic brain injury, the causal relationships between myelin disruption and the multitude of events associated with axonal damage is still unclear. Computational models of the axon can potentially clarify the damage-causality chain, provided that a mechanical description of myelin is obtained.

Here, we aim to describe membrane-component of myelinated axon at molecular level and to elucidate if and where membrane rupture can occur as a result of an injury. To achieve this, we have used molecular-based model for myelin sheath and for node of Ranvier in line with experimental composition (Table 1) and a coarse grained (CG) description. For myelin, we have built ex-novo the model based on experimental lipid composition and saturation from central nervous system (O’Brien and Rouser, 1964; Manzoli et al., 1970; Inouye and Kirschner, 1988b; Bosio et al., 1998). For the node of Ranvier, we have used an already available plasma membrane model (Ingólfsson et al., 2014), since too few experimental data on the axolemma composition for central nervous system are available. Molecular dynamics (MD) simulations have been used to characterize the molecular system at equilibrium and under deformation. Finally the deformation results from molecular simulations are combined with the axonal deformation obtained using an axonal finite element (FE) model (Montanino et al., 2020).

2 Methods

2.1 Molecular Model

2.1.1 Model for Myelin Sheath

The myelin is described as lipid bilayer. The membrane model contains 17 different types of lipids, distributed asymmetrically between two leaflets, labelled as extracellular and cytoplasmic leaflet. The extracellular leaflet is characterized by having phosphatidylcholine (PC) (7%), phosphatidylethanolamine (PE) (8%), cholesterol (CHOL) (43%), cerebroside (39%), sphingomyelin (SM) (2%), phosphatidylserine (PS) (1%) and other anionic lipids, while the cytoplasmic leaflet contains PC (11%), PE (27%), CHOL (44%), SM (3%), PS (12%) and other anionic lipids. Note that the reported values are in mol%. Supplementary Table S1 lists all the lipid types used and their disposition in the leaflets together with the level of saturation of fatty acids. The membrane model was built accounting the available experimental data: we use the values derived by Inouye and Kirschner (Inouye and Kirschner, 1988b) for myelin in central neurons system to define the lipid disposition in the extracellular and cytoplasmic leaflet; the work of Manzoli and coworkers (Manzoli et al., 1970) to define the saturation of the lipid tails for phospholipids, while the saturation values for cerebroside have obtained from O’Brien and Rouser (1964) and Bosio et al. (1998).

2.1.2 Model for Node of Ranvier

To describe the region of Node of Ranvier, we use a lipid bilayer embedded with ion channels. A channel concentration of 3,086 channels/was used to describe the node. The mammalian plasma membrane model designed by Ingólfsson et al. (2014) and deposited on MARTINI webpage (http://www.cgmartini.nl/) was used to describe the membrane. The model contains 63 different types of lipids distributed asymmetrically between two leaflets. The extracellular leaflet has a higher level of tails saturation and contains PC (36%), PE (6%), CHOL (31%), SM (19%), glycolipids (6%), and other lipids (2%). The cytoplasmic leaflet, which has a higher level of polyunsaturation, contains PC (17%), PE (25%), CHOL (29%), SM (9%), PS (11%), anionic phosphatidylinositol (PI) (2%), and other lipids (7%).

As model for the protein we used the Homo sapiens Nav1.1. The three-dimensional (3D) structure of Nav1.1 has not been yet resolved experimentally, however, the encoding gene is known (SCN1A gene) (Malo et al., 1994). We used the 3D structure obtained previously using homology modeling (for details see (Montanino et al., 2020)). As template, the cryo-electron microscopy structure of putative sodium channel from American cockroach, NavPaS (PDB ID: 5X0M) (Shen et al., 2017) was used. The structure had 100% confidence (matching probability) and 48% sequence identity (identical residues) with NavPaS. The missing terminal domains and linkers were added and the structure was energy minimized and shortly equilibrated in the PC lipid bilayers at the atomistic level using CHARMM36 force field (Pastor and MacKerell, 2011; Best et al., 2012) before building the CG model.

2.1.3 Force Field

The membrane and protein systems were described at the coarse-grained level using the MARTINI2.2 force field (Marrink et al., 2004; Marrink et al., 2007; Monticelli et al., 2008; de Jong et al., 2013) together with one-bead non-polar water model (Marrink et al., 2007). In the MARTINI model, small groups of atoms (3–4 heavy atoms) are united into beads which interact with each other by means of empirical potentials.

2.2 Systems Setup

2.2.1 Myelin Model

The myelin was modelled as a bilayer having the composition reported in Table 1. The initial lipid coordinates were constructed using the 2,000 lipids and an initial box of 24 × 24 × 12 nm. The equilibrated 20 nm bilayer (at 16 μs) was used as template to build a 40 nm bilayer and a multi-bilayers system. The 40 nm bilayer system has a total of 8,000 lipids and was placed in a cubic box of 41 × 41 × 14 nm and solvated by about 138,000 CG water beads. The multi-bilayers system was built using a distance of 1.7 and 1.6 nm between extracellular layers and cytoplasmic layers, respectively, and was placed in a cubic box of 21 × 21 × 48 nm and solvated by about 106,000 CG water beads. To all the systems NaCl was added to mimic the ionic strength at physiological condition (150 mM NaCl).

2.2.2 Node-of-Ranvier Model

A channel concentration of 3,086 channels/μm2 was used to describe the node. The bilayer (containing a total of 6,500 lipids) was placed in a cubic box (42 × 42 × 18 nm) and solvated by about 190,000 CG water beads. NaCl was added to mimic the ionic strength at physiological condition (150 mM NaCl). To achieve a concentration of 3,086 channels/μm2, 16 copies of proteins were embedded in a larger bilayer (containing 15,000 lipids in a cubic box of 69 × 69 × 20 nm). The proteins were located at 5 nm distance from each other.

2.3 Molecular Dynamics Simulations

All MD simulations were performed using the GROMACS simulation package, version 2016 (Abraham et al., 2015) (manual.gromacs.org/2016). The bilayer systems were equilibrated at constant temperature (37°C) and pressure (1 bar). The temperature was held constant using velocity rescale thermostat (Bussi et al., 2009) with a time constant of 1.0 ps. The pressure was held constant using semi-isotropic Parrinello-Rahman barostat (Parrinello and Rahman, 1981) with a time constant of 12 ps (compressibility of 3 × 10−4 bar−1). The Verlet cutoff scheme (Páll and Hess, 2013) and a timestep of 10 fs (for myelin system) and 15 fs (for protein-membrane system) were used. Periodic boundary conditions were applied. Non-bonded interactions were calculated between all beads using a cutoff of 1.1 nm. Long-range electrostatic interactions were treated using a reaction field potential (Tironi et al., 1995) with switching distance of 1.1 nm in line with MARTINI setting. After a short equilibration, 16 μs were performed for both node-of-Ranvier and 20 nm myelin models, while 25 μs MD simulations were performed for 40 nm and multi-layer myelin model. Supplementary Figure S1 shows equilibration of box dimension and energy within the first 16 μs. Note, the protein positions were kept fixed in only the simulations at equilibrium to allow the relaxation of the lipid distribution.

To evaluate membrane behavior under mechanical stress we performed 6 μs simulations at constant areal strains (NPAT ensemble) for the node-of-Ranvier and 40-nm myelin model. No position constrains was applied to the proteins Values < 0.05 were considered for areal strain. The last 4 μs are used for data production. To monitor the pore formation, we first fast pre-deformed the bilayers in x direction up to 30% strain (using a deformation speed of 1 × 10−5 nm/ps), then starting from 30% strain two independent simulations were performed using slower deformation speeds: 2.5 × 10−6 nm/ps and 25 × 10−6 nm/ps. Stretching simulations were performed for the node-of-Ranvier and myelin (40 nm and multi-layer) model. No position constrains was applied to the proteins.

2.4 Simulation Analysis

To describe the structure of lipid bilayers, we calculated bilayer thickness and number of lipid contacts. Bilayer thickness was calculated using the method proposed by Pandit et al. (2003), Pandit et al. (2004) and implemented in APL@Voro software (Lukat et al., 2013). Using Voronoi tessellation, the methods identifies the normal distance between phosphorus atoms from two leaflets that are “vertical neighbors” of each other.

Lipid contacts were calculated by counting the number of neighboring lipids within 1.5 nm of proteins or lipids by considering the first tail bead after the headgroup of lipid or polar group of cholesterol molecules. Enrichment factors were calculated as (Corradi et al., 2018):

In which and where x was chosen 0.7 and 2.1 nm for lipids around the protein in node-of-Ranvier model and 1.5 nm for lipid distribution around each lipid type in myelin model.

The partition coefficient was calculated as /, where [solute] denotes the concentration of water molecules in the lipid bilayers and water solution, respectively. was 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.

The software VMD (Humphrey et al., 1996) was used for graphical representations. Reported values were averaged on 5 μs and the errors were obtained by dividing the data production into five parts and calculating the standard error between them, if it was not specified differently.

2.5 Mechanical Properties

Area compressibility modulus () is defined as the derivative of surface tension as function of the areal strain,where γ is the surface tension and the areal strain () is defined as . Eq. 2 was chosen to be in line with experimental approach (Needham and Nunn, 1990; Rawicz et al., 2000) used on biological membrane. Evans et al. (1976) and Needham and Nunn (1990) observed a linear relation between surface tension for small areal change (<0.05 areal change) in their study on red blood cell and chlosterol-rich membrane, respectively.

To calculate area compressibility modulus, we used NPAT simulations at three strain values ( < 0.05). We divided production data in 4 parts and calculated the average γ and values γ is calculated according to the followingwhere , , and are diagonal elements of pressure matrix and is the box height along z.

We performed linear regression between the obtained γ and values, the slope of the regression line is and the standard error of the slope is the standard error of .

2.6 Multiscale Approach

To provide insights into the initiation of axonal damage, we combined the FE model of a generic portion of an axon with the molecular-based myelin and node-of-Ranvier models. The node-of-Ranvier model was used to exemplify the approach in Figure 1. Firstly, the axon FE model was utilized to simulate typical stretch injury scenarios, then as a result of axonal deformation, maximum local deformations happening at the cortex level was extracted and applied to the molecular models. The multiscale approach together with details on axon FE simulations was previously described and published in Montanino et al. (2020). The axon FE model (Figure 1C) was validated and described in Montanino and Kleiven (2018). The axon FE model is a 8 μm long representative volume of an axon consisting of three main compartments: a microtubule (MT) bundle, the neurofilament network and the axolemma-cortex complex wrapping the entire structure. To create a general axonal behavior, 10 different FE axon models were generated by randomly moving the MTs discontinuities locations, while keeping the average MTs length of the original model. Details regarding the element and material modelling choices can be found in Montanino and Kleiven (2018). All simulations were performed in LSDYNA using an implicit dynamic solver and the 1st and 2nd principal strains (, ) in the cortex plane were extracted as function of axonal strain for axonal strain rates of 1, 20, and 40/s.

FIGURE 1

Representation of the combined finite element-molecular based modeling approach (Montanino et al., 2020). (A) Finite element axonal model where a quarter of the model was removed to reveal the arrangement of the inner structures. (B) Single undeformed cortex element and a juxtaposed virtual lipid bilayer. (C) molecular bilayer model. (D) The axon finite element model is deformed of a quantity (). (E) Maximum localized deformations () are extracted and applied to the molecular model. (F) Visualization of the induced pore upon deformation of the molecular model.

3 Results and Discussion

3.1 Structural Features of Myelin Model

The myelin membrane was modelled as a lipid bilayers, having 17 different types of lipids, distributed asymmetrically between two leaflets, labelled as extracellular and cytoplasmic leaflet (Figure 2). The lipids distribution and saturation were assigned accounting for a collection of experimental data on leaflet mole composition, fatty acid length and saturation. In particular, we used the values derived by Inouye and Kirschner (1988b) for myelin in central neurons system to define the lipid disposition in the extracellular and cytoplasmic leaflet; the work of Manzoli et al. (1970) to define the saturation of the lipid tails for phospholipids, while the saturation values for cerebroside were obtained from O’Brien and Rouser (1964) and Bosio et al. (1998). Supplementary Table S1 reports details on the lipid composition of the myelin model together with the experimental reference. In Table 1 we compare the model’s lipids abundance with a wide collection of data on lipid composition observed for myelin sheath in the central nervous system of diverse species (human, bovine, rat). The myelin model has a composition in line with experimental data, small deviations are observed for PE, PS and glycolipid, whose weight fraction is between myelin and axolemma experimental composition.

FIGURE 2

Graphical representation of node-of-Ranvier bilayer model (70 nm) at 16 μs (A, B) and myelin bilayer model (40 nm) at 25 μs (C, D) together with a draw of a neuron. Both up and down leaflets are shown. Molecules are colored as follows: PC lipids are in blue, PE in yellow, CHOL in green, SM in orange, PS in purple, glycolipids in red, anionic lipids in pink, the rest of lipids in gray, protein in black VDW representation, and water is not visualized for clarity.

The myelin bilayers were simulated for 16 μs. Box dimensions and total energy converged in the first microseconds (Supplementary Figure S1). To check the convergence of the lipid distribution, we compare the numbers of lipid neighboring for the 20 and 40 nm bilayers: the values over the last 5 μs agrees within the fluctuations (Supplementary Figure S2).

The equilibrated myelin bilayer has a thickness of 4.51 nm (Table 2), larger than axolemma bilayers and comparable with the value derived using cable theory and the mean dielectric constant of squid axon myelin (around 4.5 nm) (Min et al., 2009). Table 3 reports the lipid enrichment factor. In the extracellular leaflet a depletion of phospholipids (in particular PC and PE, and SM) is observed around glycolipids. A deplection of PI around PI and enrichment of PS around SM are also reported. In the cytoplasmic leaflet no relevant enrichment/depletion is observed, most of the values are with 5% from an homogeneous distribution.

TABLE 2

Bilayers structural and mechanical properties. In parenthesis standard error.

Properties Simulated Models	Area compressibility Modulus K_A(mN/m)	Thickness Coefficient (nm)	Water partition logK_mem/wat	Rupture Strain
model for unmyel.axon	352 (3) a	4.168 (0.001) a	−1.75 (0.01) a	0.47 a ; 0.51 a
node-of-Ranvier model	252 (20)	4.032 (0.001)	−1.50 (0.01)	0.36 b ; 0.40 c
myelin model	339 (12)	4.506 (0.015)	−1.98 (0.01)	0.70 b ; 0.82 c

Lowest observed rupture strains are reported for each deformation rate.

Data from Montanino et al. (2020).

deformation rate 2.5 × 10−6 nm/ps.

deformation rate 25 × 10−6 nm/ps.

TABLE 3

Lipd enrichment factor in myelin bilayer. Number of lipids within 1.5 nm of a reference lipid type are reported for the extracellular and cytoplasmic leaflet separately. Standard errors are less than 0.1

Lipids type	PC	PE	PS	SM	PI	Glycolipid	CHOL
Extracellular leaflet
PC	1.15	1.16	1.15	1.11	1.14	0.84	1.08
PE	1.16	1.14	1.08	1.03	1.07	0.85	1.07
PS	1.12	1.04	0.89	1.17	0.91	0.88	1.07
SM	1.13	1.05	1.23	1.08	1.13	0.88	1.06
PI	1.13	1.06	0.94	1.10	0.38	0.89	1.07
Glycolipid	0.77	0.78	0.83	0.79	0.82	1.03	1.07
Cytosplasmic leaflet
PC	0.99	1.00	1.03	1.00	1.00	-	1.00
PE	1.00	0.99	1.04	0.97	1.04	-	0.99
PS	1.02	1.02	0.95	1.02	0.94	-	1.00
SM	1.03	0.99	1.06	1.01	1.02	-	0.98
PI	1.01	1.04	0.96	1.01	0.92	-	0.99

3.2 Node-of-Ranvier Model

Node of Ranvier was modelled as lipid bilayer with embedded ion channels. We use a plasma membrane lipid bilayers composition. Table 1 compares the composition with lipid abundance observed for axolemma in mammalian neuron system. The values are closer to the composition of axolemma than to the one of myelin. The larger deviation from experimental values is observed for the glycolipids. We have also to note that we can not exclude that the experimental composition for axolemma may be contaminated by myelin lipids (DeVries et al., 1981). In comparison with the composition used to describe myelin sheath, model for node-of-Ranvier is characterized by a lower content of glycolipids and CHOL molecules and higher content of PC and SM lipids, in line with the trend shown by experimental values. The model is also characterized by having shorter long-chain fatty acids than myelin model in line with experimental observation (DeVries et al., 1981).

To describe the proteins in the node, we use sodium channel protein type subunit alpha (Nav1.1), a characteristic protein of nodal portions of myelinates axons, belonging to the family of voltage-gated ion channels (Duflocq et al., 2008). It is known that several different ion channels are present, but currently information are available to model solely sodium channel protein type subunit alpha (Nav1.1). 16 proteins were embedded in the bilayers (equivalent to 3,000 channels/μm2) to mimic a value of ion channel concentration in line with the experimental observation for the node of Ranvier (Wann, 1993; Hu and Jonas, 2014). Below we compare the results with previously performed simulations of a bilayer, having the same lipid composition as the node but only one embedded protein, representing a condition closer to axolemma in unmyelinated axon (Montanino et al., 2020).

We checked the convergence of lipid distribution around the proteins by monitoring the lipid-protein contacts (Supplementary Figure S3). In general, the lipid-protein contacts are in overall agreement with the contact recorded for one protein embedded in the bilayer. The node-of-Ranvier bilayer has a thickness of around 4.0 nm (Table 2). High protein density makes the bilayer slightly thinner: a values of 4.2 was calculated for the bilayer with one embedded protein. The difference in lipid composition (both in tail length and head-group) between the node-of-Ranvier and myelin bilayers results difference in the thickness and the water permeability, given the proportionality between the partition coefficient and the permeability (). The myelin is thicker and is less permeable to water than the node of Ranvier. Now we can proceed to evaluate what is the effect of the composition on the mechanical feature.

3.3 Bilayer Models Under Deformation

To study the bilayers response to mechanical deformation, we performed simulations at three different areal strains. Figure 3 reports the surface tension as a function of the areal strain for each membrane model. In general, the surface tension increases linearly for small areal strain. The myelin model shows the steepest slope of surface tension-areal strain (Figure 3). That means that the myelin model is more resistant to change its area than node-of-Ranvier model.

FIGURE 3

Surface tension as function of the areal strain for myelin and node-of-Ranvier model, colored in red and black respectively. Reported data (circles) corresponds to the average on 1 μs time-window. The slope of the line is .

The area compressibility modulus, , was estimated for small areal strain values (less than 0.05) and reported in Table 2. A larger values is observed for the myelin model than for node-of-Ranvier models: a value of 339 mN/m and 252 mN/m, respectively for myelin and node-of-Ranvier bilayers, was obtained. At this point we can not distinguish if the difference is due to the protein concentration, lipid composition or a combination of them. The value for the axolemma model in unmyelinated (plasma membrane with one embedded protein) is in line with the area compressibility modulus measured for red blood cell membranes (375 ± 60 mN/m) at 37°C by Waugh and Evans (1979). Unlucky a validation for all the models toward experiments is not possible due to lack of experimental data. A previous study (Saeedimasine et al., 2019) shows that atomistic and CG simulations give different values, but similar trend upon change in lipid composition. That gives us trust in the observed trend, even if the absolute values can not be validated.

Simulations at different deformation rate were performed to identify the rupture points. We define the rupture point when at least one pore is formed in the bilayer. To detect the formation we monitored the surface tension: a jump in the surface tension corresponds to the formation of at least one pore. In the myelin model pores were observed at 0.7–0.8 strain while in the node-of-Ranvier model at around 0.4 (Figure 4). Poration in both bilayers occurs in region lacking protein and glycolipids. This type of lipids are mainly found in the outer leaflet of membranes and are known to make the bilayer more resistant to deformation thanks to interactions between sugar headgroups (Saeedimasine et al., 2019).

FIGURE 4

Node-of-Ranvier (A) and myelin (B) bilayers at 0.36 and 0.70 strain, respectively. For color code see Figure 2.

To provide insights into the initiation of axonal damage, we combined the FE model of a generic portion of an axon with the molecular-based myelin and node-of-Ranvier models. Firstly, the axon FE model was utilized to simulate typical stretch injury scenarios, then as a result of axonal deformation, maximum local deformations happening at the cortex level was extracted and applied to the molecular models (Figure 1). Figure 5 shows the relation between the applied axon deformation and the resulting maximum local deformation observed at different strain rates. We found that higher strain rates do not lead to higher local maximum strains for axonal strains less than 12%. This is due to the viscoelastic properties of tau proteins and of the cortex itself resulting in less cortex deformation at strain rate 40/s compared to 1 or 20/s (Montanino et al., 2020).

FIGURE 5

Axonal deformation versus maximum local strains at three strain rates: 1 s−1 (black), 10 s−1 (red) and 40 s−1 (green). Values obtained from finite element simulations of axonal stretch (Montanino et al., 2020). The dotted line corresponds to rupture strain observed for myelin model (and dashed-dot lines to the myelin multi-layer model), while the dashed line corresponds the to rupture strain observed for the node-of-Ranvier model.

We connect the bilayer strain at which local events occurs in molecular model (e.i pore formation) with the maximum local strain observed in the finite element simulations during axonal deformation. This allows to understand which axonal strain correspond to the bilayer rupture. The node-of-Ranvier model can withstand an applied deformation up to 36%, corresponding to axonal deformation of 9–11% (see dashed line in Figure 5). At higher strains pore formation is observed. While poration occurs at a larger cortex strain (70%) for the myelin model and this corresponds to an axonal deformation of 15–23%.

Myelin sheath consists of not one bilayer but of repeating units of double bilayers separated by aqueous layers. To verify the effect of the multi layers structure on the rupture threshold, we have built a double bilayers myelin model (Figure 6), equilibrated and deformed it as it was done for the single bilayer model. The distances between the bilayers were taken from X-ray diffraction data (Caspar and Kirschner, 1971; Kirschner et al., 1989): 1.6 and 1.7 nm were reported as widths of the spaces between membranes at the cytoplasmic and extracellular appositions, respectively, for the central nervous system. Starting from a strain value of 0.66 pore formation can occur in one of the bilayers (Figure 6). The first pore formation is followed by formation of multiple pores in different layers (at strain 0.73). This corresponds to an axonal strain of 14–22%, in line of what observed for the single-bilayer myelin model.

FIGURE 6

Graphical representation of multi-layer myelin at 25 μs. The multi-bilayers system was built using a distance of around 1.7 nm between extracellular and cytoplasmic layers. For lipid molecules color code see Figure 2 and for bilayer composition see Table 2.

All in all, our results show that mechanoporation occurs at 49–74% higher strain for the myelin and at 15–23% lower strain at the node of Ranvier compared with axolemma with one embedded protein. In vivo and ex vivo experiments on uniaxially oriented neuronal cell and spinal nerve showed injurious changes above 10% applied strain (Singh et al., 2015; Nakadate et al., 2017). Experiments conducted on non-oriented neuronal culture indicated that mechanoporation occurs only when strains higher than 10% were uniaxially applied to the culture (Hemphill et al., 2011). If mechanico-poration trigger axonal brain injury, our results indicate that such events in myelined axon will occur at the level of the node of Ranvier. Indeed, according to our models, poration in the node-of-Ranvier occurs at axonal strains going from 9% to 11%, while myelin rupture occurs at axonal strains larger than 15%.

Although our modeling approach brings some insights into the axonal injury mechanism, note that we do not aim to describe the injury event, since those events occur at time scale (fraction of seconds) not accessible to current molecular simulations. Moreover, the simulated bilayers are simplified models of cellular membrane: they do not account for all the different embedded proteins, they are described by a coarse grained model (where a group of atoms is described by a particle). Finally, the model does not account for possible protein structural change at different strains, since the protein is described as a semi-rigid body.

4 Conclusion

In this study, we combine CG molecular dynamics simulations and FE method to improve our understanding of what can trigger axonal brain injury. To achieve this we used molecular models to describe sub-cellular components of the myelinated axon, in particular myelin sheath and node of Ranvier, and a FE simulations of the axon.

The ex-novo built myelin model well reproduces the myelin experimental composition and structural feature of myelin. The model is more glycolipid-rich with longer fatty acid tails than the model used to describe the node-of-Ranvier. The myelin bilayer has a higher thickness, rupture point and area compressibility modulus than node of Ranvier. The results support a scenario where node-of-ranvier is most vulnerable component, myelin has more of a supporting role given that it likely will not porate before the axolemma, in unmyelined axon, and node-of-Ranvier do.

Combining the results with finite element simulations of the axon model, we provided quantitative evidences that mechanoporation of the axon membranes is an event that cannot be excluded in a typical axonal injury scenario. Our results indicate that mechanoporation may occur at the node of Ranvier at the thresholds proposed in the literature for axonal injury, while no myelin rupture is observed at the injury threshold. Poration occurs at more than 50% higher strain for the myelin compared with node-of-Ranvier and axolemma (in unmyelinated axon). The findings contribute to increase our knowledge of axonal sub-cellular components and help to understand better the mechanism behind axonal brain injury. Finally, our work shows the power of combining FE and MD to describe complex biological scenario, that requires the description of different length scale.