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⁻⁴ bar⁻¹). 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 (NP _z AT 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⁻⁵ nm/ps), then starting from 30% strain two independent simulations were performed using slower deformation speeds: 2.5 × 10⁻⁶ nm/ps and 25 × 10⁻⁶ 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.

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): E n r i c h m e n t ( L ) = R a t i o ( L ) x R a t i o ( L ) b u l k (1)

In which Ratio ( L ) x = ( no . L ) x ( tot .no .lipids ) x and Ratio ( L ) b u l k = tot .no . ( L ) tot .no .lipids 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 K m e m / w a t was calculated as [ s o l u t e ] m e m / [ s o l u t e ] w a t , where [solute] denotes the concentration of water molecules in the lipid bilayers and water solution, respectively. [ w a t e r ] m e m 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.

Area compressibility modulus ( K A ) is defined as the derivative of surface tension as function of the areal strain, K A = ( ∂ γ ∂ ϵ A ) T (2) where γ is the surface tension and the areal strain ( ϵ A ) is defined as ϵ A = ( A A 0 ) − 1 . 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 NP _z AT simulations at three strain values ( ϵ A < 0.05). We divided production data in 4 parts and calculated the average γ and ϵ A values γ is calculated according to the following γ = l z ( P z z − P x x 2 − P y y 2 ) (3) where P x x , P y y , and P z z are diagonal elements of pressure matrix and l z is the box height along z.

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

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 ( ϵ x , ϵ y ) in the cortex plane were extracted as function of axonal strain for axonal strain rates of 1, 20, and 40/s.

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.

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.

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/μm²) 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 K m e m / w a t and the permeability ( P m ). 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.

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.

The area compressibility modulus, K A , was estimated for small areal strain values (less than 0.05) and reported in Table 2. A larger K A 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 K A 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).

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).

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.

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.