Adhesion Process of Biomimetic Myelin Membranes Triggered by Myelin Basic Protein

The myelin sheath—a multi-double-bilayer membrane wrapped around axons—is an essential part of the nervous system which enables rapid signal conduction. Damage of this complex membrane system results in demyelinating diseases such as multiple sclerosis (MS). The process in which myelin is generated in vivo is called myelination. In our study, we investigated the adhesion process of large unilamellar vesicles with a supported membrane bilayer that was coated with myelin basic protein (MBP) using time-resolved neutron reflectometry. Our aim was to mimic and to study the myelination process of membrane systems having either a lipid-composition resembling that of native myelin or that of the standard animal model for experimental autoimmune encephalomyelitis (EAE) which represents MS-like conditions. We were able to measure the kinetics of the partial formation of a double bilayer in those systems and to characterize the scattering length density profiles of the initial and final states of the membrane. The kinetics could be modeled using a random sequential adsorption simulation. By using a free energy minimization method, we were able to calculate the shape of the adhered vesicles and to determine the adhesion energy per MBP. For the native membrane the resulting adhesion energy per MBP is larger than that of the EAE modified membrane type. Our observations might help in understanding myelination and especially remyelination—a process in which damaged myelin is repaired—which is a promising candidate for treatment of the still mostly incurable demyelinating diseases such as MS.

Introduction

The biological membrane is an important component of cellular function and metabolism. Investigation of biological membrane components, characterization of their physico-chemical properties and the study of their interactions with membrane binding proteins can answer important questions that are central for human health and disease. Several human diseases are directly connected to modification of cellular membranes (Evans, 1980; Maret et al., 1983; Vignini et al., 2007). Degeneration of the well-ordered myelin sheath that exists around nerve cells in the human brain, for instance, results in nerve conduction failure and neurodegeneration (Love, 2006; Weil et al., 2016). This phenomenon is called demyelination. Myelin basic protein (MBP) stabilizes the myelin membrane multi-layer and is an integral part of the myelin sheath (Boggs, 2006). Of particular relevance is the interaction of MBP with cytoplasmic leaflets of oligodendrocytes. Those cells are assembled to each other by MBP to generate a double bilayer membrane unit which envelopes the axons as a multi-membrane stack forming compact myelin (Quarles et al., 2006). Multiple sclerosis (MS) is a chronic inflammatory disease of the central nervous system, correlated with demyelination through membrane de-adhesion, swelling, and ultimately vesiculation of the myelin sheath (Weil et al., 2016). The lipid compositions of native (healthy) and modified (diseased) myelin sheaths have been investigated in an animal model (Ohler et al., 2004). While the native lipid composition occurs in healthy cytoplasmic myelin sheaths, the diseased lipid composition has been found in animals having experimental autoimmune encephalomyelitis (EAE), which is an animal model to study diseases such as e.g., MS that are connected to demyelination (Ohler et al., 2004). Structural properties of biomimetic native-like and EAE-diseased membranes have been investigated in vitro in a biosynthetic approach as oriented membrane bilayers or in the vesicle state (Min et al., 2009, 2011; Shaharabani et al., 2016, 2018; Raasakka et al., 2017). Recently, we found a direct connection between the bending stiffness of large unilamellar vesicles (LUV) having either EAE-diseased or native-like lipid composition and the formation of multilamellar structure which is mediated by the binding strength of MBP with the respective membrane types (Krugmann et al., 2020b).

While most of the above-mentioned studies have investigated structural aspects of those biomimetic myelin-like model membranes in their equilibrium state, knowledge about molecular properties that regulate the kinetics of the assembly process of the myelin sheath is currently still largely missing. A detailed molecular understanding of the formation of the myelin membrane systems is required both from a fundamental biological and biophysical point-of-view as well as to develop novel approaches for medical treatment of neurodegenerative diseases using a process known as remyelination. Remyelination is a process during which a degenerated myelin membrane is replaced by a new sheath, which seems to be thinner than the original, but it is still capable to maintain biological function of myelin (Prineas and Connell, 1979; Franklin and Ffrench-Constant, 2008). There are high expectations to utilize this natural process for the treatment of human demyelinating disorders such as, for instance, MS. Therefore, many efforts are made to understand the molecular mechanism of this process in detail. Currently, knowledge regarding the details of the molecular process of remyelination using vesicle fusion is still missing.

In our study, we investigated molecular aspects of the myelin formation process by using biomimetic membrane systems. We first produced a supported membrane bilayer (SMB) mimicking native or diseased-like cytoplasmic myelin on a silicon substrate and coated it with a dense MBP layer on top. The adsorption kinetics of LUV having the same lipid composition as those bilayers were then studied using neutron reflectometry (NR) as a function of incubation time. NR experiments performed under steady-state conditions at different stages of the assembly process allowed us to gain structural information of the oriented membrane systems orthogonal to the membrane plane with resolution on the nm length-scale. The kinetic adhesion mechanism of the LUV could be modeled by a random sequential adsorption (RSA) process. As the maximal adsorption area of RSA process is known as well as the bending rigidity of the LUV membrane, we could calculate the interaction energy per MBP protein. This parameter was found to be larger for native lipid composition as compared to the EAE-diseased membranes. A perspective for a future application of our approach would be to form myelin-like multilayers by LUV fusion, which would be an interesting mechanism to externally trigger remyelination by LUV fusion.

Materials and Methods

Materials

Porcine brain lipids (PC, PE, PS, and SM) and ovine cholesterol were purchased from Avanti Polar Lipids (Alabaster, AL, USA) and bovine MBP from Sigma-Aldrich (St. Louis, MO, USA).

Liposome Preparation

Porcine brain lipids (PC, PE, PS, and SM) with chain compositions given on the Avanti webpage¹ and ovine cholesterol have been separately dissolved in chloroform, or bought in chloroform. The lipids are mixed in ratios found in native and EAE modified myelin (Table 1; Krugmann et al., 2020b) and afterwards the chloroform is evaporated using a gentle nitrogen stream followed by vacuum annealing at 50°C overnight. The lipids are then dissolved in D₂O-buffer [150 mM NaCl, 10 mM 3-(N-morpholino) propanesulfonic acid (MOPS)]. The solution is shaken until the lipid cake is no longer stuck at the glass surface. If necessary remaining lipid cake is removed by pipetting. In the following the solution is sonicated for 30 min at 40°C and freeze-thawed 5 consecutive times. To remove remaining large aggregates or giant vesicles the solution is centrifuge-filtered through a 0.45 μm for 10 min at 10,000 g-force. Finally, the solution is extruded through a 100 nm membrane 21 × at 50°C.

TABLE 1

Table 1. Lipid compositions of native and EAE modified cytoplasmic myelin membranes (Krugmann et al., 2020b).

Supported Bilayer Preparation

Polished silicon blocks were cleaned by subsequent immersion in chloroform, acetone, ethanol and water under sonication for 20 min each, followed by UV/ozone irradiation for 30 min and rinsing in H₂O to verify hydrophilicity. All silicon substrates were stored in water until utilization to prevent contamination. For the neutron reflectometry measurements, each substrate was sealed in a solid/liquid cell. In detail, the bottom part of the cell is equipped with a PEEK through, which represents the liquid reservoir, while the upper part is an aluminum/Teflon lid. The silicon block is sandwiched between the two parts and the trough (reservoir) filled with Milli-Q water by injection through Teflon tubes (0.5 mm inner diameter), closed by valves. Prior to liposomes solution injection, the cell was pre-heated at 50°C to promote vesicle fusion. When the liposomes were in contact with the hydrophilic SiO surface, a supported membrane bilayer (SMB) was formed and residual liposomes were removed by rinsing with buffer solution.

NR Protocol

Neutron reflectivity measurements were performed on the neutron reflectometer D17 at the Institut Laue-Langevin (ILL, Grenoble, France) (Saerbeck et al., 2018; Krugmann et al., 2020a) operated in time-of-flight mode. The instrumental resolution ΔQ/Q was varying between 2 and 10% along the full q-range (0.006 0.35 A⁻¹) with accessible wavelengths from 2 to 30 Å and two angles of incidence (0.8 and 3.2°).

The footprint length (in the beam direction) and width (in the perpendicular direction) were L = 60 mm and W = 35 mm, respectively.

The samples were deposited at the solid/liquid interface of a silicon substrate sealed inside a solid/liquid cells with 1 mL liquid reservoir which was kept at 50°C to promote vesicle fusion; during the rest of the experiment a constant temperature of 37°C was maintained via a thermostatic bath. Contrast variation was carried out via an automatized pumping system at a flow rate of 1 mL/min for a full volume exchange of 20 mL. Firstly, the silicon blocks are characterized in three contrasts (D₂O-buffer, H₂O-buffer, silicon matched water (SMW)-buffer (38% D₂O). Afterwards the respective SMB is formed using the protocol described above. Now, the SMBs are characterized in the same three contrasts. We inject MBP in the cells to coat the SMB as described in Krugmann et al. (2020b). To remove remaining MBP solution the cells are flushed. Again, the membrane is measured at the three contrasts. Finally, the vesicles are injected and the adsorption kinetics is monitored (in D₂O). To enable very fast acquisition the beam divergence is increased and a time resolution of 1 min can be achieved (Cubitt et al., 2015).

NR Data Analysis

For the specular NR, the data first has been reduced using the COSMOS tool in the LAMP software of the ILL. The fitting of those reduced datasets was done using the Anaklasis package based on python (Koutsioubas, 2021). This software allows co-refining multiple datasets as e.g., the same sample at different contrasts or the same sample in D₂O at different steps of the kinetic process. The model assumed is a stratified layer model, while an additional adjustable multiplicative parameter is used for the proper scaling of the reflectivity curves. During data reduction a normalization factor is set accounting for neutron beam attenuation as defined by a direct beam measurement. In the pristine membrane samples only a small deviation to unity is observed which could be explained by a not perfect sample alignment. In the later steps the deviation from unity increases especially in the samples with attached vesicles. This might be accounted to diffuse scattering in these samples. Therefore, less signal is scattered specular and a value smaller than unity needs to be chosen as scaling factor. For the bare membrane we assume the layer order head|chain|head|. This layer sequence will just be called membrane in the following. When adding MBP an additional layer is added and the order is changed to membrane|MBP. After injection of the vesicles the order is changed to membrane|MBP|membrane|diffuse vesicle layer. The diffuse vesicle layer accounts for the 3D structure of the attached vesicles and has a SLD of 0.7 · 10⁻⁶ Å⁻² between the one of the chain and the head section. As front end we use in all cases Si|SiO|buffer and as back end buffer. Obtained membrane thicknesses are compiled in Table 2.

TABLE 2

Table 2. Coverage of the second membrane bilayer η, thicknesses of the membrane parts, MBP concentrations in the protein layers, and MBP layer thickness d_MBP.

Random Sequential Adsorption Simulation

The RSA simulation is modified from Erban and Chapman (2007) for our system. As starting point we take a test volume (x-y-z = 1 μm – 1 μm – 1 mm) with a constant concentration of vesicles in z-direction. In practice a test area (x-y) is defined. The x-y-surface at z = 0 is the adhesive surface. Then we define vesicles in z-direction. Since the concentration should be uniform in x-y-direction we only have to calculate the diffusion of vesicles in z-direction and only initialize the x-y-position in the case that the vesicles touch the surface.

Now we start the diffusion simulation. Therefore, we calculate the diffusion of each vesicle in time steps Δt by

$\begin{array}{ll}{z}_{new}=z_{old}+\sqrt{2D\Delta t}·\chi & (1)\end{array}$

Here, D is the diffusion constant and χ is a normal distributed random number.

In the case that the vesicles diffuse to the adhesive surface and no other adhered vesicle is overlapped in x-y-position they get adsorbed by a certain probability and are removed from the simulation. This is done in the simulation in the case that one of two conditions are fulfilled. The first condition is that the end z-position is negative. If this condition is fulfilled the vesicle adsorbs with the probability:

$\begin{array}{ll}p=P·\sqrt{\Delta t}& (2)\end{array}$

In Erban and Chapman (2007) P describes a positive constant which can be related to the rate constant of the chemical reaction between the virus surface and diffusing polymers. As the geometry of our system is similar in our case P is related to the interaction of vesicles and the MBP coated membrane surface. The second condition checks if the vesicles z-position was negative during the step. The adsorption probability is given by:

$\begin{array}{ll}p=exp\left(-z_{old}\frac{z_{new}}{D\Delta t}\right)·P·\sqrt{\Delta t}& (3)\end{array}$

At the beginning only, the vesicle concentration close to the surface gets perturbed by the surface adhesion. Therefore, only vesicle diffusion close to the surface is calculated in the beginning. The affected z-range L(t) increases is given by:

$\begin{array}{ll}L(t)=2·{\text{erf}}^{-1}0.99·\sqrt{tD}& (4)\end{array}$

In this model we assumed the vesicles with diameter d = 100 ± 30 nm as spherical hard shells which exclude the area $A_{ex}={\frac{\pi }{4}d}^2$ when attached to the membrane. The simulation has been run 10 consecutive times and the results have been averaged. The time step was set to Δt = 1 s. We estimated a value $P=0.015·1/\sqrt{t}$ leading to a quite convincing agreement with the measured data. The maximum coverage of such a system has been determined before (Wang, 1994; Cieśla and Nowak, 2016; Cieśla and Ziff, 2018) to be 55.47% for monodisperse disks—which is a good description of the 2D-projection of the vesicles on the membrane surface. In case of polydispersity this value increases as vesicles with smaller diameter can fit in smaller holes (Adamczyk et al., 1997). We simulated the effect for a polydispersity of 0.3 (Supplementary Figure 1) which we measured for vesicles of the same compositions and almost identical extrusion procedure via SANS and SAXS (Krugmann et al., 2020b). The coverage seems to be almost completely stable over time—at a time range relevant for the incubation time of the vesicles until the static measurements were started (native ~20 h, EAE modified ~10 h). After 1 day we simulated a steric vesicle coverage of 0.60 ± 0.05 (the error was estimated from the incubation time range of the samples) which is used in the further calculations as explained in the following chapter. Here, we used a time step Δt = 10 s in the simulation. The concentration of vesicles was chosen to be 1 mg/ml. From this value the number concentration was calculated by estimating the weight of one vesicle.

Accessing the Adhesion Constant k_a

Free Energy of a Vesicle

Let us consider a vesicle of area A, supposed to be constant upon deformation, adsorbed onto a plane. Its free energy reduces to two terms:

$\begin{array}{ll}{E}_{free}=E_a+E_{b}with\{\begin{array}{c}{E}_a=-k_a\cdot {\alpha }_{f}A\\ E_b=\frac{\kappa }{2}{\displaystyle \int {\left(c_1\left(s\right)+c_2\left(s\right)\right)}^2}ds\end{array}& (5)\end{array}$

The first term E_a is the adhesion energy potential, with k_a the adhesion constant and α_f the area-fraction of the vesicle surface that is flat and adheres to the plane. This term decreases with increasing adhesive-area fraction α_f.

The second term E_b is the bending elastic potential that can be viewed as the variance of the local curvatures c₁(s) and c₂(s) of the surface element ds, with κ the mean-curvature modulus of the bilayer. This term increases as the shape deviates from the sphere, i.e., for increasing flat part given by α_f.

For a given set (α_f, k_a, κ), the vesicle adopts the shape with local curvatures (c₁(s), c₂(s)) allowing $\int {(c_1(s)+c_2(s))}^{2}ds$ to be minimized with the constraint of curvature-continuity at the boundary of the flat area. Since this constraint depends only on α_f, so does the least-energy shape of the vesicle. Determining this least-energy surface is the first step for the calculation of E_free. For this we followed the procedure described in Koutsioubas et al. (2017). The fundamental assumption is the following: the least- energy surface of an adsorbed vesicle is a solid of revolution with a generatrix that is itself a least-energy curve. The main advantage of this method is to reduce the problem to a 1D-problem that can be exactly solved with a minimum computing-time consumption. Actually, the mean-curvature modulus k_c can be independently measured or kept from the literature. So that, E_free can be finally computed in the plane (k_a, α_f) (see Supplementary Figure 2). For a given adhesion constant k_a, vesicles adopt the value for the adhesive-area fraction α_f allowing E_free to be minimized, i.e., in the valley of Supplementary Figure 2. Conversely, by measuring α_f one can determine the corresponding value for k_a. The bending rigidity κ of native and EAE-diseased membranes was determined by NSE previously (Krugmann et al., 2020b) and the values of κ are compiled in Table 3. Obtained values of α_f and k_a determined in this study are summarized in Table 3 as well.

TABLE 3

Table 3. Values of the bending rigidity κ, the membrane coverage η, the fraction of the vesicles which are flat α_f, the calculated adhesion energy k_a, the area per MBP molecule A_MBP and the adhesion energy per MBP molecule ϵ_MBP.

Combining RSA Simulations and Neutron Reflectivity Measurements to Access α_f

“From the top” point of view (RSA):

In the RSA model, the addition of a new particle is sensitive to the steric hindrance of particles already adsorbed. Let us denote A_ex the maximum extension area parallel to the plane of a vesicle, i.e., its projection area on the plane, and define the “steric coverage density” ρ as:

$\begin{array}{ll}\rho =N·\frac{A_{ex}}{A_{tot}}& (6)\end{array}$

where N is the number of adsorbed vesicles and A_tot the total area of the plane.

“From the bottom” point of view (neutron reflectivity):

On the other side, neutron reflectivity is sensitive to the area of vesicle that adheres and is directly in contact with the plane. Let us denote A_f this area with no curvature per vesicle and define the “adhesive coverage density” η as:

$\begin{array}{ll}\eta =N·\frac{A_f}{A_{tot}}& (7)\end{array}$

Dividing Equation (6) by Equation (7) gives:

$\begin{array}{ll}\frac{\rho }{\eta }=A_{ex}/A_f& (8)\end{array}$

By definition, A_f writes:

$\begin{array}{ll}{A}_f={\alpha }_{f}A& (9)\end{array}$

But actually, the flatter the vesicle, the wider it is, so that A_ex also depends on α_f. For a spherical vesicle of radius R, one has $A_{ex,a_f=0}=\pi R^2$ and for a completely flattened vesicle $A_{ex,a_f=0.5}=A/2=2\pi R^2=2A_{ex,{\alpha }_f=0}$. The computation of the least-energy surface as a function of α_f shows that the departure from the linearity between these two extremities is negligible (see Supplementary Figure 3).

This leads to

$\begin{array}{ll}{A}_{ex}\simeq A·(1+2{\alpha }_f)/4& (10)\end{array}$

Equations (8) and (10) give: ρ/η = (1 + 2α_f)/4α_f leading to

$\begin{array}{ll}{\alpha }_f\simeq \frac{1}{4\left(\frac{\rho }{\eta }\right)-2}& (11)\end{array}$

Considering the situation at time t → ∞ with ρ_∞ = 0.60 ± 0.05 is estimated from the simulations, one can deduce α_f from the neutron reflectivity measurement of η_∞.

Results and Discussion

To precisely analyze the interaction of MBP coated membranes with other membranes we choose to rely on neutron reflectometry (NR). Information is obtained perpendicular to the membrane plane by detection of the specular NR signal in combination with in-plane resolution by also measuring off-specular NR data (Zhou and Chen, 1995; Jablin et al., 2011). This is shown schematically in Figure 1. The reflected neutron beam is mostly scattered specular—meaning that the final angle Φ_fi is equal to the initial angle Φ_in. The specular signal is determined by the structure perpendicular to the surface and has a resolution in the order of a few Ångstrom. However, if there is also in plane structure in the sample—in most cases roughness—a part of the beam is reflected at a different angle. This fraction is called off-specular and is only sensitive to parts of the sample which have in-plane structure.

FIGURE 1

Figure 1. Model of structures seen by the specular (A,B) and off-specular (C) NR signal. In the specular signal the layer structure in z direction is measured very precisely. However, it is more sensitive for the structure close to the silicon surface. From the fit of the specular data the coverage of the second bilayer can be estimated (top right panel). In the off-specular signal the roughness of the layer is measured (C). If this signal appears one can conclude that the layer structure is not flat.

Preparation of the Membrane System and Steady-State Characterization

The system of interest is a supported membrane bilayer (SMB) coated with MBP which is kept in a liquid cell. Into this cell large unilamellar vesicles (LUV) are injected which can diffuse toward the SMB surface and adhere to it. The main purpose of our study was to investigate the mechanism of that adhesion process. Prior to that we validated the structural properties of the membrane components—an essential prerequisite for the interpretation of the kinetic NR experiments. The different stages of the preparation process are schematically depicted in Figure 2 alongside with the average thicknesses of the membrane components.

FIGURE 2

Figure 2. Schematic representation of structural models for the different stages of the experiment. The averaged thickness values of the membrane components are indicated in the figure panels. (A) Supported membrane bilayer formed by vesicle fusion on silicon wafer. (B) MBP-coated membrane formed after injection of protein solution and incubation. (C) Double-layer membrane structure connected by a dense MBP phase in between that is formed upon adhesion of the LUVs.

As a first step the SMB was created by vesicle fusion on a hydrophilic silicon substrate (Figure 2A; Stroumpoulis et al., 2006; Rondelli et al., 2017). The neutron reflectometry curves were measured on the neutron reflectometer D17 at the ILL (Saerbeck et al., 2018) at three different contrasts having buffer compositions of 100% D₂O, 100% H₂O and a mixture (38%/62% v/v) D₂O/H₂O which matches the SLD of silicon being abbreviated as silicon matched water (SMW) in the following. The NR data are shown for a native and EAE modified SMB in Figures 3A, 4A, respectively. The fit of a simple head—chain—head layer sequence can reproduce the measured NR data nicely (see Figures 3A, 4A). A SMB with native-like lipid composition and two SMBs with EAE modified lipid composition referred to as Native01, Modified01 and Modified1 SMBs were characterized. The obtained thicknesses of the head (d_head) and chain sections (d_chain) of those three SMBs are reported in Table 2 and the respective water penetration values are given in Supplementary Table 1. The respective SLD profiles are shown in Figure 3B of the Native01 membrane and in Figure 4B of the Modified01 SMB.

FIGURE 3

Figure 3. Neutron reflectivity curves and respective SLD profiles of the Native01 sample. (A,B) The pure native membrane without MBP, (C,D) the supported membrane with MBP protein layer on top that was formed by addition of protein solution with 0.1 mg/ml MBP, and the double-membrane structure with dense intramembrane MBP-layer that was formed after vesicle adhesion in (E,F). In (F) the SLD profiles of the pristine membrane and the MBP coated membrane are shown as comparison. NR data was measured on the D17 neutron reflectometer.

FIGURE 4

Figure 4. Neutron reflectivity curves and respective SLD profiles of the EAE Modified01 sample. (A,B) The pure EAE membrane without MBP, (C,D) the supported EAE membrane with MBP protein layer on top that was formed by addition of protein solution with 0.1 mg/ml MBP, and the double-EAE-membrane structure with dense intramembrane MBP-layer that was formed after vesicle adhesion in (E,F). In (F) the SLD profiles of the pristine membrane and the MBP coated membrane are shown as comparison. NR data was measured on the D17 neutron reflectometer.

As a next step MBP was injected in the liquid cell at either a concentration of c_MBP = 0.1 mg/ml for the Native01 and Modified01 SMBs or at 1 mg/ml for the Modified1 SMB. Due to electrostatic interaction, it interacts with the negatively charged membrane and forms a dense protein layer on top of the membrane (see Figure 2B). The data of the Native01 membrane with 0.1 mg/ml MBP in combination with a fit is shown in Figure 3C. The respective data of EAE modified membrane Modified01 with 0.1 mg/ml MBP including theoretical fit is shown in Figure 4C. Those datasets have been fitted as a layer sequence head—chain—head and an additional dense MBP layer on top of the membrane with thickness d_MBP,SBL and concentration c_MBP,SBL. During these fits the fitted values of d_head and d_chain, water penetration and roughness between the layers from the SMB without MBP were used and constrained in ±5% of the best value around those values, while the thickness of the MBP layer d_MBP,SBL (constrained for Modified01: 60–100 Å, otherwise: 70–100 Å) and its water penetration were fitted. The water penetration of the chain section was left free between 0 and 5%. The SLD of MBP was estimated with the primary structure of bovine MBP using the Biomolecular Scattering Length Density Calculator² considering 90% exchange of labile hydrogen atoms. These values are in D₂O ${\rho }_{D_{2}O}=3.45·10^{-6}{A}^{-2}$, in SMW ${\rho }_{SMW}=2.55{·10}^{-6}{A}^{-2}$, and in H₂O ${\rho }_{H_{2}O}=2.00·10^{-6}{A}^{-2}.$ The respective SLD profiles of the Native01 and Modified01 SMB are depicted in Figures 3D, 4D. The fit results are reported in Table 2, the corresponding buffer penetration values are given in Supplementary Table 1. We got a convincing fit with a MBP layer with thickness between d_MBP,SBL = 60–80 Å (the determined parameters are characterized by large error bars, due to the dilute nature of the MBP layer) and a protein concentration in the range of c_MBP,SBL = 3–7% vol/vol (corresponding to MBP concentrations of around 22 and 51 mg/ml, respectively, when a protein specific volume of 0.73 ml/mg is considered) for both membrane types and both injected MBP concentrations of 0.1 and 1 mg/ml. The MBP layer thicknesses and concentrations as reported in this study are in agreement with previous observations (Krugmann et al., 2020b).

Vesicles were injected onto the MBP/SMB systems and the interaction process was studied with kinetic NR acquisition. The kinetic data will be discussed in the next section, but beforehand we take a look at the steady state specular NR measurement after an almost constant state is reached after 1–2 h. A schematic visualization of that steady-state membrane system in shown in Figure 2C. Experimental NR data of native and EAE modified conditions under steady-state conditions are shown in Figures 3E, 4E including theoretical fits. We can see a strong change in the reflectivity profiles in comparison to the membranes without adhered vesicles. As a fit model we assume a double-bilayer system because the vesicles diffuse to the membrane and adhere to it. In this case, the specular NR profile should be a double bilayer as specular NR only is sensitive for the z-structure close to the Si-surface since the penetration depth of the neutrons under such small angle is quite low. However, to account for the part of unfused adsorbed vesicles that extends into the solution, we add in our model a long 800 Å diffuse vesicle layer as last layer with a SLD of 0.7 · 10⁻⁶Å⁻² and high roughness and buffer penetration. The thickness of the MBP layer in the double bilayer model has been fitted to d_MBP,DBL ≈ 30 Å which is close to what we have measured before in between 2 bilayers in small-angle X-ray scattering experiments (Krugmann et al., 2020b). The fitted values of d_head and d_chain, water penetration and roughness between the layers from the SMB without MBP were again used and constrained in ±5% of the best value around those values. The thicknesses of the second membrane layer (d_chain, d_head,inner, and d_head,outer) are fixed to the values d_chain and d_head,outer of the first layer. They are fixed to d_head,outer as they are both not attached to the silicon substrate and therefore less ordered. Fitted parameters are the buffer penetration value corresponding to the MBP concentration between the bilayer c_MBP,DBL. The fitted buffer penetration values are given in Supplementary Table 1 and the MBP concentrations in c_MBP,DBL in Table 2. The SLD-profiles we attain form those fits are shown in Figures 3F, 4F.

The fitted structural models of the three stages of the experiment are summarized visually in Figure 2 including the averaged thickness values of the membrane components. Overall, our characterization by NR of the assembly process under steady-state conditions yields a coherent picture of the assembly process. MBP that is injected as a rather dilute concentration interacts strongly with the first SMB and forms a concentrated fluid phase with concentration in the range of 3–7% and a thickness of 60–80 Å. Upon formation of the second membrane layer the MBP layer thickness is reduced to around half of its initial thickness (~30 Å) and a concentration of 0–1% vol/vol for the 0.1 mg/ml samples and 11% for 1 mg/ml sample. This MBP layer is, hence, regulating the assembly and fusion process of the LUV with the SMB. However, due to the complexity of the system in its final state and the dilute nature of the MBP layer, it appears that our modeling is sensitive only to the thickness of the MBP layer, leading to relative underestimation on protein concentration.

Furthermore, from the buffer penetration value h_chain,2 of the second chain section (given in Supplementary Table 1) it is possible to calculate the coverage η of the second bilayer as the hydrophobic chain section is nearly water free in both cases. Therefore, it is possible to calculate the coverage as η = 1 − h_chain,2. The definition of the coverage η can be expressed mathematically as the area of the second bilayer A_dbl divided by the total area A_tot of the first bilayer of the sample

$\begin{array}{ll}\eta =\frac{A_{dbl}}{A_{tot}}.& (12)\end{array}$

In Table 2 the coverage values of η are given for the investigated membrane systems.

In Figure 5A the 2D off-specular map of an EAE-modified sample in the final steady-state of LUVs bound to the MBP-coated SMB measured at the MARIA reflectometer is given. For technical details of the instrument MARIA, we refer to the instrument paper (Mattauch et al., 2018). In Figure 5B the 2D image of an EAE modified sample of LUVs bound to the MBP-coated SMB measured using the D17 reflectometer is shown. Clearly, an off-specular signal is visible in both cases. This off-specular signal is proof that the samples have strong surface roughness which is an indication for vesicles adhered to its surface (Ott and Kozhevnikov, 2011).

FIGURE 5

Figure 5. Off-specular signal of vesicle coated membranes seen by (A) classical monochromatic and (B) TOF neutron reflectometer. In (A) an exemplary measurement at the classical monochromatic instrument MARIA (Mattauch et al., 2018) from an EAE-modified sample of our previous paper (Krugmann et al., 2020b) is shown. And in (B) an exemplary measurement of an EAE modified membrane at the TOF reflectometry mode using an incident polychromatic neutron beam of the D17 instrument is shown as comparison.

Kinetics of the Vesicle Adhesion Process

Time-resolved NR experiments were carried out to understand the process of adhesion. In Figure 6A an exemplary kinetic NR measurement of a modified membrane is shown. To maximize the NR signal and to reduce the incoherent background, time-resolved NR experiments have been performed only at one contrast in 100% D₂O buffer. Clearly, in the time range of one to a few hours the adhesion process is happening. The off-specular data is growing in this time as well. The specular curves in Figure 6A can be fitted by the earlier explained double bilayer model. Here, only the coverage (η = 1 − h_chain,2), water penetration of the diffuse vesicle layer its roughness and the scaling factor are fitted, while the rest of the remaining parameter values are fixed to the ones obtained from the static measurements of the membrane with adhered vesicles. The influence of the head section is not that strong as the contrast to D₂O is higher to the chain section. In Figure 6B the respective SLD profiles are depicted. In Figure 6C the normalized coverage (scaled to the RSA simulation curve) from the specular fit, the normalized off-specular signal (scaled to the RSA simulation curve with subtracted background) and the behavior represented by the RSA simulation curve are plotted together. The important observation here is that the normalized experimental coverage values obtained from both specular and off-specular NR data are identical within the statistical uncertainty and their time-dependence can be reproduced by the theoretical model of an RSA simulation.

FIGURE 6

Figure 6. Kinetics of the reflectivity measurements: In (A) a waterfall plot of the kinetic reflectivity measurements of the EAE modified membrane with 1 mg/ml MBP is shown (sample Modified1). Experimental NR data are represented by symbols as a function of incubation time. Solid lines are theoretical fits using the double-membrane bilayer model with intramembrane-MBP layer and diffuse vesicle layer for 3D structure of the attached vesicles. Structural parameters of the membrane profile were kept fixed during the individual fits and only the coverage fraction of the second membrane layer (and the diffuse vesicle layer water content and roughness) was fitted. In (B) the respective SLD profiles are shown. In (C) the kinetics of the off-specular signal of the three membranes (empty symbols), the kinetics of the coverage fit parameter of the three membranes obtained from specular NR data (filled symbols) and an RSA simulation of vesicles with 50 nm radius and 15 nm polydispersity (solid black line) are shown. In (D) the log-log-plot of the vesicle coverage ρ(inf)−ρ(t) over the time is shown for the off-specular experimental values and the RSA simulation. The thin line indicates power-law behavior with a slope of −1/2, which is the theoretical prediction of an RSA process of flat disks on a surface (Adamczyk et al., 1997).

RSA is a model in which particles diffuse and adhere to a surface and cannot detach from it. The details of the RSA simulation are explained in the Materials and Methods section. Strikingly, the kinetics of the specular data—representing the formation of a double bilayer—and the off-specular data—representing the adhesion of vesicles to the surface—follow the same trend and can both be described by the same RSA simulation. We can conclude here that both signals are governed by the same kinetic time constant of the diffusion driven adhesion of the vesicles. The deformation of the vesicles which is responsible for the double bilayer formation is, therefore, quasi-instantaneous in comparison to the adhesion process. In RSA the vesicles that adsorb to the membrane are described as disks that—when adsorbed—exclude a certain area A_ex. After waiting for an infinite amount of time (t = inf) the coverage of vesicles ρ(t)–which is different from η as not the whole excluded diameter is in contact with the membrane—a coverage ρ(inf) is reached. To verify if RSA is the correct kinetic process to describe our data we plot the logarithm of ρ(inf)-ρ(t) over the logarithm of time in Figure 6D. This should give a slope of −1/2 for an RSA process of disks on a surface, which is indeed the case of our NR data when we exclude the first initial time points (Hinrichsen et al., 1986).

Calculation of Adhesion Energy per MBP

From the above-mentioned results, it is possible to estimate the adhesion energy per area k_a. This is achieved by minimizing the free energy which is constituted of the bending energy and the adhesion energy (see section Free Energy of a Vesicle). By accounting for some geometrical assumptions (see section Combining RSA simulations and neutron reflectivity measurements to access α_f) we can estimate the adhesion energy per area k_a of the native and modified membranes with 0.1 mg/ml MBP. Values of k_a are given in Table 3. With this we can calculate the interaction energy per MBP molecule ϵ_MBP.

$\begin{array}{ll}{ϵ}_{MBP}=k_a·A_{MBP},& (13)\end{array}$

where A_MBP is the area per MBP molecule we take from Krugmann et al. (2020b) calculated with the MBP molecule volume from Stadler et al. (2014). Input values of A_MBP and the obtained values of ϵ_MBP are compiled in Table 3 for the native and EAE-diseased membranes. We can see that the so-estimated adhesion energy per MBP-molecule in native membranes—calculated by this model—is larger as in EAE modified membranes which might explain the formation of a stable myelin sheath. The A_MBP values calculated from fits of D17 data that are presented in this publication have large statistical uncertainty, but they are in agreement with values of A_MBP that we have reported previously (Krugmann et al., 2020b). The large error bars of the D17 fit parameters A_MBP (Table 3) lead to high errors of ϵ_MBP as well. Therefore, it is not possible to draw any conclusions on the difference between native and EAE-diseased membranes using these values.

In the case of the EAE modified membrane with 1 mg/ml MBP concentration, we have a coverage fraction η_dbl of around 39%. This would lead to a high value of k_a in the range of the value of silicon $k_a(Si)=0.5-1mJ/m^2$ where we already observe vesicle fusion triggered by the strong adhesion (Anderson et al., 2009). Therefore, we believe that it is unphysical and probably we have partial vesicle fusion occurring. This case is not covered in our model and would of course increase η_dbl without the vesicles needing to be flattened as much. Therefore, the calculated k_a would decrease dramatically. Partial vesicle fusion might also happen in the lower protein concentration case where the values are close to those measured in silica. In that case our calculated values of k_a and ϵ_MBP might be overestimated.

Conclusions

In this manuscript, we report on the adhesion mechanisms of LUV onto a supported membrane that was coated with an MBP layer. Structural characterization at the different stages of membrane assembly process has been performed under steady-state conditions, while the vesicle fusion process has been monitored as function of the incubation time. The overall aim was to investigate differences between native and EAE modified biomimetic cytoplasmic myelin membranes and to gain information on the myelin sheath assembly process at a molecular scale. From the steady-state experiments we demonstrate that a concentrated MBP layer is formed on top of a single SMB with a protein concentration of around 3–7%. Adhesion of LUV with local formation of a second membrane bilayer results in an MBP layer that is reduced in thickness. Concerning the LUV adhesion kinetics, we show that the adhesion of the vesicles can be modeled by an RSA simulation. In addition, we provide evidence that the vesicles are reshaping after the adhesion. Their contact surface with the membrane is maximized until the elasticity of the membrane leads to a force in opposite direction and an equilibrium state is established. Our measurements show that this process depends on both membrane type and MBP concentration, while its kinetics is entirely limited by the diffusion of the vesicles toward the membrane surface agreeing with the observed RSA mechanism. In combination with the bending rigidity of those membranes, which we investigated in earlier studies (Krugmann et al., 2020b), we can provide an estimate of the interaction force per MBP molecule: In the native-like membrane case it is larger than in the EAE-diseased membranes. Our new results can also help with the interpretation of our last paper (Krugmann et al., 2020b). In that study, we have seen that LUV having EAE modified lipid composition tend to stick to each other by a strong attractive interaction being mediated by the binding strength of MBP. As we show in here, this seems not to be because of a stronger interaction of MBP with these EAE-diseased membranes as the adhesion energy per MBP-molecule is higher for native membrane than for the EAE-diseased membrane. Instead, the driving force of that effect appears to be due to the higher MBP concentration in the MBP layer on top of EAE modified membranes as compared to that of native-like membranes.

We found indications that partial vesicle fusion is occurring when EAE membranes were treated with the higher MBP concentration of 1 mg/ml. By applying external stimuli like osmotic shock or temperature shocks it might be possible to achieve complete vesicle fusion resulting in a perfect double bilayer system with an intercalated MBP layer in between the two bilayers. This would be a significantly improved model system to study the initial growth process of the myelin sheath from a molecular perspective, and it would allow us to study basic biophysical properties of artificial biomimetic nerve cells in vitro.

Data Availability Statement

The datasets analyzed for this study can be found in the ILL data repository (data.ill.eu) (Krugmann et al., 2020a).

Author Contributions

AS, AR, AK, and SF designed the research. BK, AK, LH, and SM performed the research. BK and AK analyzed the data. DL contributed analytical modeling. BK and AS wrote the manuscript with input from all other authors. All authors read and approved the finalized manuscript.

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.

Acknowledgments

This work was based upon experiments performed on the instruments MARIA operated by Jülich Centre for Neutron Science at the Heinz Maier-Leibnitz Zentrum, Garching, Germany and on D17 at the Institut Laue-Langevin, Grenoble, France. The authors acknowledge the PSCM at the ILL for the use of the preparation laboratories and for the preliminary sample characterization during the beamtime.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2021.631277/full#supplementary-material

Footnotes

1. ^https://avantilipids.com/ (accessed March 1, 2021).

2. ^Biomolecular Scattering Length Density Calculator. http://psldc.isis.rl.ac.uk/Psldc (accessed March 17, 2021).

References

Adamczyk, Z., Siwek, B., Zembala, M., and Weroński, P. (1997). Influence of polydispersity on random sequential adsorption of spherical particles. J. Colloid Interface Sci. 185, 236–244. doi: 10.1006/jcis.1996.4540

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, T. H., Min, Y., Weirich, K. L., Zeng, H., Fygenson, D., and Israelachvili, J. N. (2009). Formation of supported bilayers on silica substrates. Langmuir 25, 6997–7005. doi: 10.1021/la900181c

PubMed Abstract | CrossRef Full Text | Google Scholar

Boggs, J. M. (2006). Myelin basic protein: a multifunctional protein. Cell. Mol. Life Sci. CMLS 63, 1945–1961. doi: 10.1007/s00018-006-6094-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Cieśla, M., and Nowak, A. (2016). Managing numerical errors in random sequential adsorption. Surf. Sci. 651, 182–186. doi: 10.1016/j.susc.2016.04.014

CrossRef Full Text | Google Scholar

Cieśla, M., and Ziff, R. M. (2018). Boundary conditions in random sequential adsorption. J. Stat. Mech. Theory Exp. 2018:043302. doi: 10.1088/1742-5468/aab685

CrossRef Full Text | Google Scholar

Cubitt, R., Saerbeck, T., Campbell, R. A., Barker, R., and Gutfreund, P. (2015). An improved algorithm for reducing reflectometry data involving divergent beams or non-flat samples. J. Appl. Crystallogr. 48, 2006–2011. doi: 10.1107/S1600576715019500

PubMed Abstract | CrossRef Full Text | Google Scholar

Erban, R., and Chapman, S. J. (2007). Time scale of random sequential adsorption. Phys. Rev. E 75:41116. doi: 10.1103/PhysRevE.75.041116

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans, W. H. (1980). Hepatic plasma-membrane modifications in disease. Clin. Sci. 58, 439–444. doi: 10.1042/cs0580439

CrossRef Full Text | Google Scholar

Franklin, R. J. M., and Ffrench-Constant, C. (2008). Remyelination in the CNS: from biology to therapy. Nat. Rev. Neurosci. 9, 839–855. doi: 10.1038/nrn2480

PubMed Abstract | CrossRef Full Text | Google Scholar

Hinrichsen, E. L., Feder, J., and Jøssang, T. (1986). Geometry of random sequential adsorption. J. Stat. Phys. 44, 793–827. doi: 10.1007/BF01011908

CrossRef Full Text | Google Scholar

Jablin, M. S., Zhernenkov, M., Toperverg, B. P., Dubey, M., Smith, H. L., Vidyasagar, A., et al. (2011). In-plane correlations in a polymer-supported lipid membrane measured by off-specular neutron scattering. Phys. Rev. Lett. 106:138101. doi: 10.1103/PhysRevLett.106.138101

PubMed Abstract | CrossRef Full Text | Google Scholar

Koutsioubas, A. (2021). Anaklasis: A Compact Software Package for Model Based Analysis of Specular Neutron and X-Ray Reflectometry Datasets. Journal of Applied Crystallography. Available online at: https://github.com/alexandros-koutsioumpas/anaklasis (accessed April 1, 2021).

Koutsioubas, A., Appavou, M. S., and Lairez, D. (2017). Time-resolved neutron reflectivity during supported membrane formation by vesicle fusion. Langmuir 33, 10598–10605. doi: 10.1021/acs.langmuir.7b02459

PubMed Abstract | CrossRef Full Text | Google Scholar

Krugmann, B., Gutfreund, P., Haris, L., Koutsioubas, A., Micciulla, S., Radulescu, A., et al. (2020a). Kinetics of Myelin Membrane Fusion Triggered by Myelin Basic Protein. Grenoble: Institut Laue-Langevin (ILL). doi: 10.5291/ILL-DATA.8-02-870

PubMed Abstract | CrossRef Full Text

Krugmann, B., Radulescu, A., Appavou, M. S., Koutsioubas, A., Stingaciu, L. R., Dulle, M., et al. (2020b). Membrane stiffness and myelin basic protein binding strength as molecular origin of multiple sclerosis. Sci. Rep. 10:16691. doi: 10.1038/s41598-020-73671-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, S. (2006). Demyelinating diseases. J. Clin. Pathol. 59, 1151–1159. doi: 10.1136/jcp.2005.031195

CrossRef Full Text | Google Scholar

Maret, A., Potier, M., Salvayre, R., and Douste-Blazy, L. (1983). Modification of subunit interaction in membrane-bound acid β-glucosidase from Gaucher disease. FEBS Lett. 160, 93–97. doi: 10.1016/0014-5793(83)80943-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Mattauch, S., Koutsioubas, A., Rücker, U., Korolkov, D., Fracassi, V., Daemen, J., et al. (2018). The high-intensity reflectometer of the Jülich Centre for Neutron Science: MARIA. J. Appl. Crystallogr. 51, 646–654. doi: 10.1107/S1600576718006994

PubMed Abstract | CrossRef Full Text | Google Scholar

Min, Y., Alig, T. F., Lee, D. W., Boggs, J. M., Israelachvili, J. N., and Zasadzinski, J. A. (2011). Critical and off-critical miscibility transitions in model extracellular and cytoplasmic myelin lipid monolayers. Biophys. J. 100, 1490–1498. doi: 10.1016/j.bpj.2011.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Min, Y., Kristiansen, K., Boggs, J. M., Husted, C., Zasadzinski, J. A., and Israelachvili, J. (2009). Interaction forces and adhesion of supported myelin lipid bilayers modulated by myelin basic protein. Proc. Natl. Acad. Sci. U.S.A. 106, 3154–3159. doi: 10.1073/pnas.0813110106

PubMed Abstract | CrossRef Full Text | Google Scholar

Ohler, B., Graf, K., Bragg, R., Lemons, T., Coe, R., Genain, C., et al. (2004). Role of lipid interactions in autoimmune demyelination. Biochim. Biophys. Acta 1688, 10–17. doi: 10.1016/j.bbadis.2003.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Ott, F., and Kozhevnikov, S. (2011). Off-specular data representations in neutron reflectivity. J. Appl. Crystallogr. 44, 359–369. doi: 10.1107/S0021889811002858

CrossRef Full Text | Google Scholar

Prineas, J. W., and Connell, F. (1979). Remyelination in multiple sclerosis. Ann. Neurol. 5, 22–31. doi: 10.1002/ana.410050105

CrossRef Full Text | Google Scholar

Quarles, R. H., Macklin, W. B., and Morell, P. (2006). “Myelin formation, structure and biochemistry,” in Basic Neurochemistry, Molecular, Cellular and Medical Aspects, 7th Edn., ed S. Brady (Amsterdam: Academic Press), 51–71.

Google Scholar

Raasakka, A., Ruskamo, S., Kowal, J., Barker, R., Baumann, A., Martel, A., et al. (2017). Membrane association landscape of myelin basic protein portrays formation of the myelin major dense line. Sci. Rep. 7:4974. doi: 10.1038/s41598-017-05364-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Rondelli, V., Brocca, P., Tranquilli, N., Fragneto, G., Del Favero, E., and Cantù, L. (2017). Building a biomimetic membrane for neutron reflectivity investigation: complexity, asymmetry and contrast. Biophys. Chem. 229, 135–141. doi: 10.1016/j.bpc.2017.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Saerbeck, T., Cubitt, R., Wildes, A., Manzin, G., Andersen, K. H., and Gutfreund, P. (2018). Recent upgrades of the neutron reflectometer D17 at ILL. J. Appl. Crystallogr. 51, 249–256. doi: 10.1107/S160057671800239X

CrossRef Full Text | Google Scholar

Shaharabani, R., Ram-On, M., Avinery, R., Aharoni, R., Arnon, R., Talmon, Y., et al. (2016). Structural transition in myelin membrane as initiator of multiple sclerosis. J. Am. Chem. Soc. 138, 12159–12165. doi: 10.1021/jacs.6b04826

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaharabani, R., Ram-On, M., Talmon, Y., and Beck, R. (2018). Pathological transitions in myelin membranes driven by environmental and multiple sclerosis conditions. Proc. Natl. Acad. Sci. U.S.A. 115, 11156–11161. doi: 10.1073/pnas.1804275115

PubMed Abstract | CrossRef Full Text | Google Scholar

Stadler, A. M., Stingaciu, L., Radulescu, A., Holderer, O., Monkenbusch, M., Biehl, R., et al. (2014). Internal nanosecond dynamics in the intrinsically disordered myelin basic protein. J. Am. Chem. Soc. 136, 6987–6994. doi: 10.1021/ja502343b

PubMed Abstract | CrossRef Full Text | Google Scholar

Stroumpoulis, D., Parra, A., and Tirrell, M. (2006). A kinetic study of vesicle fusion on silicon dioxide surfaces by ellipsometry. AIChE J. 52, 2931–2937. doi: 10.1002/aic.10914

CrossRef Full Text | Google Scholar

Vignini, A., Nanetti, L., Moroni, C., Tanase, L., Bartolini, M., Luzzi, S., et al. (2007). Modifications of platelet from Alzheimer disease patients: a possible relation between membrane properties and NO metabolites. Neurobiol. Aging 28, 987–994. doi: 10.1016/j.neurobiolaging.2006.05.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J.-S. (1994). A fast algorithm for random sequential adsorption of discs. Int. J. Mod. Phys. C 05, 707–715. doi: 10.1142/s0129183194000817

CrossRef Full Text | Google Scholar

Weil, M.-T., Möbius, W., Winkler, A., Ruhwedel, T., Wrzos, C., Romanelli, E., et al. (2016). Loss of myelin basic protein function triggers myelin breakdown in models of demyelinating diseases. Cell Rep. 16, 314–322. doi: 10.1016/j.celrep.2016.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, X. L., and Chen, S. H. (1995). Theoretical foundation of X-ray and neutron reflectometry. Phys. Rep. 257, 223–348. doi: 10.1016/0370-1573(94)00110-O

CrossRef Full Text | Google Scholar