MEDUSA: A cloud-based tool for the analysis of X-ray diffuse scattering to obtain the bending modulus from oriented membrane stacks

An important mechanical property of cells is their membrane bending modulus, κ. Here, we introduce MEDUSA (MEmbrane DiffUse Scattering Analysis), a cloud-based analysis tool to determine the bending modulus, κ, from the analysis of X-ray diffuse scattering. MEDUSA uses GPU (graphics processing unit) accelerated hardware and a parallelized algorithm to run the calculations efficiently in a few seconds. MEDUSA’s graphical user interface allows the user to upload 2-dimensional data collected from different sources, perform background subtraction and distortion corrections, select regions of interest, run the fitting procedure and output the fitted parameters, the membranes’ bending modulus κ, and compressional modulus B.


Introduction
Cellular functions, such as mobility, division and vesicle trafficking, are intrinsically related to a cell's ability to comply to deformation [1][2][3].Especially the cell membrane's endurance against bending forces is critical for a cell's survival.The Helfrich Hamiltonian [4] for the energy of a bent symmetric membrane is: where uðrÞ describes local spatial deviation of the bilayer center in the out-of-plane direction (with respect to the membrane) and the integral is over the area covered by the membrane.The bending modulus κ is a material property.As such, it varies with temperature and the molecular composition of the membrane and is a particularly appropriate measure of overall membrane elasticity.
In XDS experiments, the sample consists of stacks of solid supported multilamellar membranes and the X-ray scattering is measured.Such samples are widely used to probe the structure of synthetic lipid bilayers [20][21][22][23][24] as well as more complex biological membranes [25,26].Importantly, such samples are ideal to analyze the mechanical properties of these molecular structures as it was first introduced by Lyatskaya et.al [12] and has been subsequently applied to single [13,14,27] and multi-component [28][29][30] lipid bilayers culminating in studies on native red blood cell membranes [15,16,31].Advancements in comprehending the molecular-level mechanics within these membranes have given rise to a variety of sophisticated techniques, facilitating the utilization of endogenous membranes for various biotechnological applications [32,33].
An exemplary result is shown in Fig 1A .The most intense scattering is specular (q k =0 Å −1 ) with peaks that originate from the average lamellar repeat distance in the stack of membranes, and it includes the sharp line of reflectivity from the solid support.The space between neighboring membranes is occupied by water and the thickness d w of this water layer is controlled by the environments relative humidity (RH).In experiments, d w is commonly calculated as d w = d − d m , where d represents the lamellar repeat spacing and d m is the membrane thickness (also known as Luzzati thickness) as defined in [34].d w is small (d w � 13 Å for RBC membranes [25]) below 98% RH and out-of-plane membrane fluctuations are suppressed.Well hydrated samples (>99.9%RH), on the other hand, have a significantly larger (d w � 28 Å) water layer between the membranes [15], enabling out-of-plane fluctuations of the individual membranes.This results in an additional diffuse off-specular scattering (see Fig 1B).Membrane fluctuations are thermally driven and their amplitude and spectrum vary with the bending rigidity κ.The intensity in the off-specular regime is closely related to the membrane structure factor as will be discussed in depth below, so one can experimentally determine κ by fitting the model structure factor as a function of κ to the experimental out-of-plane intensity profile.
Calculating the structure factor is computationally challenging.Here, we introduce MEDUSA (MEmbrane DiffUse Scattering Analysis), a cloud-based analysis tool to determine the bending modulus, κ, from XDS experiments.MEDUSA uses GPU accelerated hardware and a parallelized algorithm to run the calculations efficiently.MEDUSA's graphical user interface allows the user to upload 2-dimensional data collected from different sources, perform background subtraction and distortion corrections, select regions of interest, run the fitting procedure and output the fitted parameters, the membranes' bending modulus κ and interaction modulus B. The tool is available at https://medusa.genapp.rocks.
As diffuse scattering is typically orders of magnitude weaker than specular scattering, these measurements require X-ray sources with sufficient intensities.Both synchrotron sources and optimized rotating anode in-house machines have been used.A critical element of the experimental design is the humidity controlled sample environment as measurements need to be performed close to 100% RH.Various designs have been developed in the past.We will describe two successful chamber designs and provide insight into the challenges and caveats.

Membrane scattering theory
The intensity measured in an X-ray diffraction experiment can be written as where F(q z ) is the form factor and SðqÞ is the structure factor.In the case of solid supported lipid bilayer, this structure factor has the form [12]: where J 0 is the zero order Bessel function and δu n (r) is the height-height pair correlation function.The functions H z (z) and H r (r) account for experimental limitations that reduces the effective size of coherently scattering membrane patches.A limiting factor in the out-of-plane A 2-dimensional X-ray intensity map recorded on a stack of POPC bilayers.Measurements were performed 100% RH.The most intense scattering is specular (q k = 0) with additional diffuse off-specular scattering that originates from out-of-plane membrane fluctuations.B Intensity profile of A along the out-of-plane scattering vector q k at q z =3q 1 and 3.5q 1 , where q 1 refers to the 1 st order lamellar peak.C 2-dimensional map of Eq 4. D Simulated line-cuts of C at q z =3q 1 and 3.5q 1 . https://doi.org/10.1371/journal.pcbi.1011749.g001 direction is the number of membranes that can be stacked with a sufficient degree of order.The in-plane direction is limited by the sample's dimensions, the footprint of the beam (typically in the order of <200 μm) on the sample, but most importantly, by the diameter of coherently scattering membrane patches.These patches are different from the lipid domains that occur in mixed lipid bilayers and need to be understood as regions that scatter coherently and are finite due both to finite beam coherence and to sample mosaicity.The membrane stack then consists of many of these regions.The structure factor for such a patchy sample has been described in detail in [12] and [35].Briefly, they are assumed to be cylindrical patches with a Gaussian distributed diameter L r (average patch size L r and variance σ r ) and a Gaussian distributed height L z (average patch size L z and variance σ z ) and define size effect functions [12]: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi XDS is the off-specular scattering observed in the diffraction experiment.Since F(q z ) only depends on q z , by Eq 2, diffuse scattering is solely governed by the structure factor S(q z , q r ) which depends on the inter-and intra-lamellar height-height pair correlation function δu n (r) of the membrane stack.This is shown in Fig 1C which depicts a 2-dimensional map of Eq 4 and illustrates the out-of plane contribution (Fig 1D ).A non-vanishing δu n (r) arises from thermally excited out-of-plane fluctuations; an analytical expression has been derived from the stack's free energy [12].

Stack model
Eq 1 describes the energy of a single lipid bilayer.For a sample consisting of a stack of bilayers, an interaction between bilayers must be added [12,36]: where B is a modulus that accounts for the interaction between neighboring membranes in the harmonic approximations.[37].
Calculating the height-height pair correlation function δu n (r) from Eq (6) has been described in detail in [37].Briefly, membrane fluctuations are governed by thermal energy and can be separated into normal modes by transforming the out-of-plane displacement u n ðrÞ into Fourier space (u n ðrÞ !U n ð QÞ).Q spans the Fourier space of the membrane fluctuations and differs form the scattering vector q.The free energy functional in Eq (6) decouples in Fourier space.The equipartition theorem then assigns 1  2 k B T of energy to each normal mode.This allows calculating the power spectrum of the membrane fluctuations.
In Fourier space, the height-height pair correlation function is proportional to this power spectrum and an analytical expression of du n ðrÞ was derived [37], where J 0 is the zero order bessel function, q 1 = 2π/d, and ξ and η are known as Caille ´parameters which relate to the bending modulus κ and the membrane interaction modulus B through [12] Z Here, k B is the Boltzmann constant and T is the temperature.One can determine both, the bending modulus κ and the interaction modulus B, independently by fitting the structure factor S(q z , q k ) to experimental data.

Data processing
MEDUSA is built using the GenApp framework [38].GenApp enables a convenient web deployment of code in any programming language.Details on GenApp can be found elsewhere [38] and the following paragraphs solely focus on the two back-end programs.

Reduction
MEDUSA supports 2-dimensional intensity maps in the Tag Image File Format (.tiff) file-format.In a first step, data are loaded and reduced by a reduction library developed in Python [39].The program flow of the data reduction is visualized in Fig 2A .Intensities can be stored on a linear or logarithmic scale.A rudimentary background-subtraction routine is implemented as discussed below.However, it is recommended to subtract any instrumental background for optimal results.Subsequently, the meridional and azimuthal angle of the diffraction experiment will be referred to as θ and X.Most 2-dimensional flat detectors subtend the spherical coordinate system of the reciprocal space spanned by θ and X and consequently measure a distorted image (see Fig 3).This distortion can be corrected when taking the geometry of the X-ray instrument and the data acquisition by the detector into account.However, MEDUSA allows the selective enabling and disabling this distortion correction to accommodate for specialized detectors.
Let us first discuss pixels located at q k = 0 Å −1 .The azimuth angle θ k for every pixel is determined by where PX size is the size of a single pixel, PX k is the k-th pixel and PX 0 is the location of the direct beam on the detector.In the same way, The out-of-plane and in-plane component of the scattering vector are subsequently calculated using ð11Þ once both angles are determined for a given dataset.The 2-dimensional X-ray data are reduced in two different ways.1.Data are cropped to a 12 pixel wide rectangular box centered at q k = 0 and averaged along q k .This provides an outof-plane, meridional, intensity scan that enables users to determine the repeat d-spacing of the stack.2. Data are cropped into several rectangular boxes each centered at different values of q z that are specified by the user.The box width is also specified by the user.One option then averages along q z within the box.This determines what will be named one q z in-plane line cut.A minimum of two such line cuts is required for determination of K C and B. Another option provides horizontal q z line cuts for each vertical pixel in the box.Finally, the 2-dimensional dataset as well as the out-of-plane and in-plane line cuts are visualized using the plotly graphing library.

Fitting
The intensities for each q z line cut are simultaneously fitted to the membrane's structure factor in Eq 2. The parameters in the fit are K C and B and a separate multiplicative factor for each line cut to take into account the q z dependence of the form factor F(q z ).Calculating S(q z , q k ) is computationally challenging and the subprogram for this was thus written in C++.The algorithm was based on a previous program [12] but was modified to allow for GPU acceleration with the Compute Unified Device Architecture (CUDA) provided by the Nvidia Corporation https://doi.org/10.1371/journal.pcbi.1011749.g003[40].GPU acceleration generally works by splitting the computation workload of a given problem between multiple processors.
The CUDA toolkit allows splitting of a processing job into threads that are grouped in blocks.The number of threads per block is a hardware specific quantity.The maximum number of blocks is independent of the hardware and is only limited by the CUDA toolkit [41].As a result, a processing job can be split into as many threads as required.The effective speed gain is limited by the hardware.A single Geforce GTX-1080 TI graphics card, for instance, has 3584 physical CUDA cores and allows 1024 threads per block.
The flow diagram of the implemented algorithm is depicted in Figs 2B and 4. The program uses the program_options toolbox from the boost C+ library to handle user input and can operate in two modes: It can calculate the 2-dimensional structure factor for a given q z and q k range, or it can fit a provided data set.Both routines rely on an algorithm that calculates the structure factor for given values q z , q k , ξ, η and q 1 .Calculating and fitting the structure factor in Eq (4) is non-trivial due to the nested integration and summation and requires computational approximations.The term that solely depends on q z can be isolated from the structure factor allowing us to rewrite Eq (4) We can consequently calculate Λ(r) only once for a given q z before solving the integration in Eq (13) numerically for the desired values of q k (hereafter referred to as q k -profile).The functions H r (r) and H z (nD) were introduced to account for the finite size of membrane domains.This is convenient as it reduces the required range for n and r [37].
The first step in calculating Λ(r) is to compute the height-height pair correlation function δu n (r).It is computationally useful to use Eq (7) for n < 30 and r < 1000 Å and employ the approximation [12] du n ðrÞ ¼ for all other values for n and r.Both equations for δu n (r) are independent of the scattering vector q.We can introduce the transformation r ¼ rx, where r is the radius for ξ = 1 and η=1.
Any other combination of (ξ, η) can then be calculated by simply rescaling r: This enables us to calculate an array of δu(r, 1, 1) at the beginning of the algorithm for logarithmically spaced floating point values 10 −4 Å< r < 10 6 Å and linearly spaced integer values 0 < n < 1000.Eq (15) is then applied to calculate an array for δu(r, ξ, η).The integration in Eq (7) is performed numerically using the adaptive.gsl_integration_qagiu algorithm provided by the GNU scientific library [42].The trapezoidal rule can not be used due to the apparent singularity in the integrand in Eq (7) (x !0).
In the same way, arrays for H r (r) and H z (nd) (see Eqs ( 4) and ( 5)) are pre-calculated.Again, logarithmically spaced floating-point values 10 −4 Å< r < 10 6 Å (10,000 values in total) and linearly spaced integer values 0 < n < 1000 were used in the calculation.The calculation of H r (r) ) before calculating a single q k profile for given values of q z , q k , ξ, η and q 1 .https://doi.org/10.1371/journal.pcbi.1011749.g004 is further accelerated using the CUDA toolkit by splitting the process into 10 blocks with 1024 threads each.Each thread then calculates the integration in Eq (5) for fixed values r and n and stores the results in an array.
In the next step, the algorithm calculates Λ(r) (Eq (12)) for integer values of −1000 < n < 1000 using the array entries from all predetermined functions.This process is split into 2 blocks with 1024 threads each.Each thread solely calculates the summation in Eq (12) for a given value r and stores the results in an array.Finally, the program calculates Eq (13).The numerical integration is performed using the trapezoidal rule with 1 Å< r < 10 6 Å and a step size of 1 Å.Values of Λ between the grid points of the predetermined arrays are determined from cubic interpolation.This process is once again parallelized.Two blocks with 1024 threads each are defined, where each thread is instructed to compute the integration for a fixed value q r .Eq (4) represents the structure factor for a finite membrane stack, but does not account for characteristics of the X-ray instrument.In a real-world experiment, the structure factor is convoluted with the beam's footprint on the sample.The beam profile in the described setup is circular with a Gaussian distribution with a a standard deviation of σ q specified by the user in both spatial directions.The determined q k profile is thus convoluted with a Gaussian distribution to account for this beam geometry.Multiple q k -profiles are calculated by looping through multiple q z to calculate a 2-dimensional scan of the structure factor S(q z , q k ).
The bending modulus κ and the membrane interaction modulus B can be determined independently from XDS data by fitting the q k dependence of the calculated structure factor S(q z , q k ) at more than one independent values of q z ¼ q ð1Þ z and q z ¼ q ð2Þ z to the experimental data.For this purpose, a Levenberg-Marquardt least squares fit was implemented using the gsl_mul-timin_fminimizer from the GNU Scientific library [42].
The function to be minimized is given by the sum of the squared residuals where Y l [k] is the interpolated value of S(q z , q k ) at discrete values q ðlÞ z and q ðkÞ r and y l [k] are the corresponding experimental values.σ l [k] are the experimental errors.

Workflow
User authentication.MEDUSA requires a user authentication via the GenApp interface [38].The login can be accessed through the link in the upper right corner and requires the user-id and the password.New users can register by clicking on the avatar icon and entering a user-id, a password (minimum 10 characters) and an email address.This authentication requirement is motivated by two considerations: First, it reduces the risk of unauthorized access and protects computational resources.Second, the data structure of GenApp differs when authentication is enabled.GenApp creates a user-specific folder to execute and store all related data within a module.This folder is set to be constant for all modules.This allows a convenient communication between the reduction and fitting module (see below).Additionally, it enables users to store fit results and revisit previous fits.

Experimental setups
Experiments can be performed using either a synchrotron source or an in-house rotating anode machine.Data should be recorded far enough in the q z direction to obtain the lamellar repeat spacing d and to see robust diffusion scattering.A typical range is for q z from 0 to 0.5 Å −1 .Data should also be recorded far enough on both sides of the meridian to include the full decay of the diffuse scattering in the q k direction in order to allow adequate background subtraction.A typical range is for q k from -0.5 to +0.5 Å −1 .Two setups have been routinely used.The instrument is equipped with a 9 kW CuKα rotating anode tube and a RIGAKU HyPix-3000 2-dimensional semiconductor detector.Multilayer optics consisting of a focusing mirror, a 5 degree soller collimator, and a 5 mm monocapillary collimator provide a circular beam with a diameter of �200 μm, a divergence of 0.008 rad and an intensity of 10 8 counts/mm 2� s.The wavelength is λ = 1.5418Å with a spread of Dl l ¼ 1%.The detector has an array of 775×385 pixels, each measuring 100×100 μm 2 .Each pixel is a single photon counter with a bitdepth of 32 bit.A beam-block was installed to attenuate the intensity from the direct, i.e. nonscattered, beam.
2) The experiments have been performed on a synchrotron source and on standard rotating anode sources using a different setup.The incident angle is controlled by rocking the sample rather than changing the angle of the incoming synchrotron beam, and the detector has been placed at a fixed position in space.Instead of taking multiple exposures for every sample angle, every pixel of the detector accumulates counts during the entire sample rotation.Importantly, while convenient in improving the signal-to noise ratio, that type of rocking leads to an additional convolution to the data analysis to account for the pixel intensity originating from different values of q.Although the MEDUSA software is not written specifically for that setup, when applied to data obtained that way, MEDUSA returns values of the moduli that agree quite well with what was obtained with software [37] specifically written for that setup.
For either X-ray optical configuration a critical component in the experimental setup is a sample chamber capable of achieving between 99.9% and 100% RH in order to fully hydrate the stack of membranes from the vapor.This requires a chamber design that minimizes condensation and ensures a homogeneous temperature within the chamber.Three such chambers have been constructed [14,40,43].
Here, we describe the chamber specifically designed for the first set-up above ( [40]) that shares features of the earlier chambers.It was machined in aluminum which is opaque for Xrays, so windows are machined on either side of the chamber and are sealed with a 13 μm thick Kapton foil.This polymeric material was chosen for its high transmittance for X-rays Table 1.A summary and description of the input parameter.

2D X-ray Data (.tiff)
Single file containing 2-dimensional scatting data in the .tiffformat q-cut parameter Specification out-plane scattering vector for each q k line-cut.
Pixel to be averaged Specification of the box width for averaging along q z in every line-cut.

Instrument Specific Parameters
Sample to Detector Distance Distance between sample and detector in mm.
Required to calculate the accurate q-space of the data-set Detector Pixel Size Dimension of individual pixels (CCD Detector or single photon counter) in mm.
Required to calculate the accurate q-space of the data-set Origin in X and Y Position of the direct beam on the detector measured in pixel.
λ Wavelength of the instrument Distortion in q par or q z Enable to apply a distortion correction in q k and q z (half-width at half-maximum �0.05 Å −1 ) which can be avoided, if desired, by using mylar instead of capton.The sample is placed on a stage in the center of the chamber.A basin at the bottom of the chamber below the sample is filled with a hydration solution.The humidity inside the chamber is controlled through the choice of salt and the salinity of the saline solution.For many synthetic membranes, ultra-pure water is advantageous as it allows reaching close to 100% RH.However, biological membranes such as RBC membranes are highly hygroscopic and can swell until they are washed off the silicon wafer, so it can be advantageous to lower the relative humidity by using a 40 mg/ml K 2 SO 4 solution for hydration.An alternative synchrotron sample design suggested in [14] was constructed at NIH in Bethesda.Fig 7C shows a schematic and Fig 7D is a photograph of the chamber.The nearly cubic aluminum chamber (5 × 5 × 6 inches) was designed for ease of use in demanding synchrotron experiments.Water circulates through channels in the one-inch thick walls, top, and bottom of the chamber from a temperature-controlled bath.The chamber has double windows with circulating helium in between, using thin mylar material to minimize scattering interference.The inner chamber contains water to hydrate the sample vapor and ensure thermal contact.Hydration is aided by a sponge on the chamber's top, increasing evaporation surface area, and a Peltier cooler cooling the sample relative to the vapor.Sample holders for flat and cylindrical substrates rotate with programmable motors.The chamber is flushed with helium to minimize air scattering and sealed off.A thermocouple monitors temperature inside the chamber.
The discussed theory is not limited to X-rays but generally applies to scattering of other particles within a wavelength range between 1 and 10 Å. Especially neutrons are of interest.Humidity chambers for neutron beam-lines have been constructed out of a single piece of aluminum (by J. Katsaras) without the need for windows which make full hydration difficult.Chambers exist in most facilities.Neutron scattering in the form of Neutron Spin Echo (NSE) spectrometry has been widely used to measure the membrane's bending rigidity.However, the established analysis has recently been challenged in complex lipid mixtures [44,45].If there is sufficient neutron intensity to observe diffuse scattering, this might provide novel insights into the membrane's mechanical properties.Importantly, samples with multiplex lipid mixtures may be selectively deuterated enabling a domain-sensitive measurement of the bending rigidity.

Availability and future directions
We have develped MEDUSA (MEmbrane DiffUse Scattering Analysis), a cloud-based analysis tool to determine the bending modulus, κ, from X-ray diffuse scattering experiments.MEDUSA uses GPU accelerated hardware and a parallelized algorithm to run the calculations efficiently.MEDUSA's graphical user interface allows the user to upload 2-dimensional data collected from different sources, perform background subtraction and distortion corrections, select regions of interest, run the fitting procedure and output the fitted parameters, the membranes' bending modulus κ and interaction modulus B. The tool is available at https://medusa.genapp.rocks.

Fig 1 .
Fig 1.A 2-dimensional X-ray intensity map recorded on a stack of POPC bilayers.Measurements were performed 100% RH.The most intense scattering is specular (q k = 0) with additional diffuse off-specular scattering that originates from out-of-plane membrane fluctuations.B Intensity profile of A along the out-of-plane scattering vector q k at q z =3q 1 and 3.5q 1 , where q 1 refers to the 1 st order lamellar peak.C 2-dimensional map of Eq 4. D Simulated line-cuts of C at q z =3q 1 and 3.5q 1 .

Fig 2 .
Fig 2. A) Flow diagram of the reduction routine (written in Python).B Flow diagram of the fitting program (written in C++ and CUDA).Critical subroutines are highlighted in light blue and light green, and are visualized in greater detail in Fig 4.https://doi.org/10.1371/journal.pcbi.1011749.g002

Fig 3 .
Fig 3. Most 2-dimensional flat detectors subtend the spherical coordinate system of the reciprocal space spanned by the angle θ.

Fig 4 .
Fig 4. Flow diagram of the subroutines.The program pre-calculates arrays of the height-height pair correlation function δu n (r) and the finite size effect functions (Eqs 4 and 5) before calculating a single q k profile for given values of q z , q k , ξ, η and q 1 .

Fig 5 .
Fig 5. Screenshots of the reduction tab GUI: A Reduction-specific and Instrument-specific parameters can be entered by the user.B The progress of the reduction is visualized in a text-field and progress-bar.Once complete, the program visualizes the 2-dimensional intensity maps, the reflectivity profile as well as the q k line cuts (source: https://medusa.genapp.rocks/medusa/).https://doi.org/10.1371/journal.pcbi.1011749.g005

Fitting ParametersξFig 6 .
Fig 6.Screenshots of the fitting tab GUI: A Fitting-specific parameters can be entered by the user.B The progress of the fitting is visualized in a textfield reporting the log-output from the GPU-accelerated program.The determined bending rigidity K c , the interaction modulus B as well as χ 2 are updated for every iteration.Once complete, the program visualizes the q k line cuts with the respective best-fit line and provides a download link to the reduced and fitted data (source: https://medusa.genapp.rocks/medusa/).https://doi.org/10.1371/journal.pcbi.1011749.g006

Fig 7 .
Fig 7. The setup of the X-ray diffraction machine is schematically sketched.The central components: X-ray tube, collimator optics, humidity chamber and detector are marked in the graphics.https://doi.org/10.1371/journal.pcbi.1011749.g007