Lipid bilayer model

[aglavic]

PR Scope

This PR deals with three additions to the model language definition and code:

1. Add the possibility to describe SubStackType classes directly in the stack string using the syntax ClassName{p1: <p1>, p2: <p2>}. This description is optional and useful for classes with a few non-default arguments.
2. A related description for named objects that are resolved later. In the future these object should have a separate database in SLDdb to be more flexible: ClassName{name} is used in the stack parameter.
3. Add a Leaflet and Bilayer class as high level description and parameterization of lipids and bi-layers used in soft-matter experiments. Both support a direct name resolution for currently hardcoded lipids that, in the future, should be migrated to a name resolution similar to materials.

The functionality is based on discussions with @maxskoda and @StephenCLHall and expanded using input from thie PR discussion.

Usage

Here an example of how the new functionality could be used to describe a sample with 3 bilayers:

minimal:

stack: "Si | 2 (Bilayer{DMPC}) in D2O | (Bilayer{DMPC}) in H2O | H2O"

full:

stack: "Si | lbl1 | lbl2 | (lbl3) in H2O | H2O"
sub_stacks:
  lbl1:
    sub_stack_class: Bilayer
    outer: PC
    inner: DM
  lbl2:
    sub_stack_class: Bilayer
    outer: PC
    inner: DM
    apm: {magnitude: 0.8, unit: nm^2}
    outer_hydration: 0.3
    inner_hydration: 0.3
    inner_hydration_2: 0.4
    outer_hydration_2: 0.5
  lbl3:
    sub_stack_class: Bilayer
    outer: 
        formula: C10H18O8NP
        number_density: 0.3135 # 1/nm^2
    inner:
        formula: C26H54
        number_density: 0.128 # 1/nm^2
    apm: 0.9 # nm^2
    outer_hydration: 0.3
    inner_hydration: 0.3
    inner_hydration_2: 0.4
    outer_hydration_2: 0.5
  materials:
      PC:
          formula: C10H18O8NP
          number_density: 3.135 # 1/nm^2
      DM:
          formula: C26H54
          mass_density: 0.779 # g/cm^3
globals:
    default_solvent:  D2O

This is an example of a solid-liquid experiment with three lipid bi-layers. In the simple example they are defined given the known DMPC lipid, The full description give explicit outer and intter materials (which could also be names to search in SLDdb). In the current implementation, the hydration referrs to the "default_solvent" but placing the item in a SubStack allows to overwrite that environment (here defined by "(...) in H2O").

It may be advantagous to add a heuristic for definition of default_solvent from the stacks first or last material, if it is a liquid. Not sure how to check for this at the moment.

The PR incldues two new examples for the full functionality (10+11) that include lipid leaflets and bilayers with various degrees of modification.

The model building blocks can be used by software to implement the model with known components. An example of that for refnx can be seen in resolve_and_plot_model.py where the Leaflet and Bilayer objects are resolved to LipidLeaflet objects instead of slabs.

Lipid resolution

Currently, the class will include a set of pre-defined Lipids. Later it is envisaged to add the ability of composite materials into the SLDdb to search for Lipids specifically.
https://github.com/refnx/refnx/blob/main/tools/app/lipid_properties.ipynb
https://github.com/reflectometry/molgroups/blob/main/molgroups/lipids.py

Guest insertion

This is another topic to address in the future. I'd leave it out of the PR, which is already pretty heavy, and address later once I could discuss it with you in more detail.

[andyfaff]

Some comments:

• unfortunately a bilayer class isn't able to handle a lipid monolayer (such that is found at the air-liquid interface). Is it better to have a class describing a leaflet, like refnx's LipidLeaflet, that can handle either case?
• It doesn't look like the class can handle asymmetric bilayers, which is an important use case. This is another argument for a LipidLeaflet-like functionality.
• I think the class should be able to handle guest insertion, either into the head or tail regions, or both. A guest molecule could be something like cholesterol, protein, antibiotic, protein, etc. The volume fraction of guest should be fittable.

[aglavic]

Hi @andyfaff ,

this is still the draft PR, I didn't finish implementing the initial class and need some more input from Max about the Math behind it. I thought having the PR already created allows the discussion to start before finalizing the coding.

How we describe these systems should be discussed between you guys with soft-matter background. I don't see an issue having both a LipidLeaflet and a Bilayer class, if that is what the community wants to use.

Insertion of guests would probably done by some optional parameter, we didn't get to that discussion as we wanted to start with the simplest case.

[maxskoda]

Dear @andyfaff and @aglavic ,

Thanks for the comments and starting on the implementation.
All @andyfaff 's comments are valid of course. I would say the following:

• the objective for me is usually ease of use, rather than completeness or versatility
• this means that I don't see it as a problem for several similar constructs to co-exist, as long as they are documented. So having LipidLeaflet should not preclude us from having "symmetricBilayer" and "asymmetricBilayer" too.
• I think the idea here was more to allow constructs other than those which resolve to layers in orsopy, i.e. to allow "components" with the minimum number of arguments (for ease of use), which then then resolve in orsopy to "component X" with some payload output parameters, which can then be used in a well-defined way by the fitting software to assemble a model, which is human readable and can be edited later by the user.

Steve and I have finally put together the following:

ORSOLipidBilayer_Example.py

This represents the call signature e.g. "symmetricBilayer((POPC, 0.8), (POPS,0.2))". Our suggestion is that:

• orsopy should recognise "symmetricBilayer" and parse its parameters (lipid types and respective volume fractions - exact call is not important, it could be two lists for lipid names and volume fractions, or a series of dicts, whatever is easiest)
• orsopy should collect the following information from the SLD database for each lipid name (if exists):
  • tailVolumes, tailScatteringLengths, headVolumes, headScatteringLengths
• orsopy should then return:
  • the above four lists,
  • a list for the volume fractions from the original input
  • additional parameters required to build the model (areaPerMolecule, headgroupHydration, bilayerCoverage, bilayerRoughness)

In the fitting software, assembling the model becomes the quite simple, provided a "symmetricBilayerORSO" class exists with the functionality listed in the above attached file:

def customLipidBilayerORSO(params, bulk_in, bulk_out, contrast):

# Parameters populated from ORSO Model Language
subRough           = params[0]
SiO2Thickness      = params[1]
SiO2Roughness      = params[2]
SiO2Hydration      = params[3]
areaPerMolecule    = params[4]
headgroupHydration = params[5]
bilayerCoverage    = params[6]
bilayerRoughness   = params[7]

oxide_slab = [SiO2Thickness, 3.46e-6, SiO2Roughness, SiO2Hydration, 1] # SiO2 SLD populated from ORSO SLD Database

lipidBilayer = symmetricBilayerORSO(tailVolumes           = [944, 944],
                                    tailScatteringLengths = [-267.3e-6, -267.3e-6],
                                    headVolumes           = [319, 312],
                                    headScatteringLengths = [602.7e-6, 845e-6],
                                    volumeFractions       = [0.8, 0.2],
                                    areaPerMolecule       = areaPerMolecule,
                                    headgroupHydration    = headgroupHydration,
                                    bilayerCoverage       = bilayerCoverage,
                                    bilayerRoughness      = bilayerRoughness,
                                    bulkSolventSLD        = bulk_out):


output = [oxide_slab] + lipidBilayer.stacks

return output, subRough

This function (RasCAL stye here) the returns the complete stack of regular layers plut bilayer layers based on the fitting parameters listed at the top.

Steve has also include a potential example, like Artur's suggestion for defining a lipid, which is not in the SLDDB. This should also work, but I am concerned that the defining the syntax correctly without an orsopy "model builder" will become trickier for the average user.

Please let us know what you think.

Thanks,

Max and Steve

[aglavic]

@maxskoda
I will try to implement that and include the new option of volume fractional mixes of lipids. For orsopy I would prefer to create new materials from voume fractions (fomula+density) and not just SLDs so it is radiation type independent.

I could include an example for parsing the model class in a way as you describe, similar to what I'm doing in genx. That should generally be indpendent of the layer generation in orsopy but lead to the same resulting sld profile.

[aglavic]

I'm now happy with the status to go for review. Please add others you think should have a look at the implementation.

I've seen the physics being encoded quite diffently in refnx then what I got from Max, so I guess there are various ways to describe the same thing. The documentation will have to clearly define the resolution of physics to model slabs, so I'd suggest a sketch like that with the associated math. Here is how I implemented and understood the various parameters with volume $V$ and hydration $h$:

$V_{i}$ = $\frac{1}{\text{number density}_{i}}$
$VF_{i}$ = $1- h_i$ (Lipid volume fraction)
$V_{total}$ = $\frac{ V_{lipid} }{VF_{lipid}} $ (Actually calculated separately for heads and tails)
$d_{heads}$ = $ \frac{ V_{heads} }{ \text{APM} \cdot VF_{heads}}$
$d_{tails}$ = $ \frac{ V_{tails} }{ \text{APM} \cdot VF_{tails}}$
$ \text{SLD}_{head-slab}$ = $ (1-h_{heads}) \text{SLD}_{heads} + h_{heads}\text{SLD}_{solvent}$
$ \text{SLD}_{tail-slab}$ = $ (1-h_{tails}) \text{SLD}_{tails} + h_{heads}\text{SLD}_{solvent}$

[bmaranville]

This seems like effectively a plugin system/contract, where orsopy recognizes certain inputs for a substack as corresponding to materials which can be looked up from an sld table, and then the results fed into a custom class indicated by the "class" attribute. Is this contract documented? It seems like in the example given symmetricBilayer the contract is expected to be hard-coded on both ends (the schema is fixed). If this is to be extensible, then both orsopy and the community of consumers have to agree on how each new class is going to be interpreted (which parts will be pre-processed by orsopy, and which be passed through directly to the consumer class).

Just as a thought, this might be more flexible if everything were passed directly to the consumer class as-is, and the consumer is free to use orsopy to look up SLD for materials as it needs to (any consumer program already requires orsopy to be installed and imported at the time of use)

[aglavic]

@bmaranville It goes into this direction, but the idea was more based on the concept we have for the ORSO header in general. The orsopy library is used to, at the same time, to define the schema of the sample model as well as an example for the implementation. What is special here is the physical interpretation of the results, which is also shipped within the orsopy classes.

Yes, this will require much more documentation to be added. Right now the PR is focused on the definition as well as implementation.

That said, a program that wants to use orsopy to analyze the model (recommanded) has four options:

1. Just load the data and interprete the information itself from the class attributes. This means, all name resolution needs to be done by the client, which can become quite complex due to the felxibility we wanted for the model.
2. Call SampleModel.resolve_stack(), this makes resolution of the stack string with names into associated class object but no flattening of SubStack objects is performed and no materials are resolved. Then the client can do all interpretation of the data a bit more easily, but sill needs to understand more about the model language, itself. In this process, some normalization of objects is also perfomed (e.g. defaults for units are resolved so that the client can always assume the attributes to be "Value" objects, not potential floats with default unit.
3. Call SampleModel.resolve_to_blocks() does all of the above but also flattens SubStacks. It keeps SubStackType classes (like the Bilayer) in tact. To allow the client to efficiently use those, these do material resolution. There might be a way to do this even later, but this is as much complexity as I could take for now considering all the material mixing necessary. In any case, the resolution of materials is coverned by MaterialResolver derived classes, that can be registered by the client to perform such analysis. By default the SLDdb based resolver is used.
4. Call SampleModel.resolve_to_layers() flattens everything to slabs. This should be the fall-back that allows any software to represent the model with individual slabs that have SLD, thickness and roughness. For software that can deal with complex objects, this fall-back also exists for the SubStackType classes, so it can be used for anything not desribed within them.

Have a look at examples/resolve_and_plot_model.py, showing how this can be done for refnx.

Edit:
Just to clearify, the thing posted by @maxskoda are not part of the orsopy library. I think it was supposed to describe the function used by a specific reflectivity software to build a model from a ORSO sample model description. We could include some exmples for other programs like rascal, refnx was just simples for me to put in the library for resolving and plotting the sample model example files.

[aglavic]

To illustrate how this works, this is the resulting model for refnx from a simple stack by running:
python resolve_and_plot_model.py ""Si | (Leaflet{rDMPC}) in D2O | Bilayer{DMPC} | H2O""

####################### Neutron Model ##################
from refnx.reflect import SLD, LipidLeaflet, ReflectModel, Slab, Structure

structure = Structure()

m = SLD((2.07370547828562-2.375752014432938e-05j))
structure |= Slab(0.0, m, 3.0, name='Si')

structure |= LipidLeaflet(
            apm=70.0,
            b_heads=(0.0006007780000000002-2.2512013099632003e-08j),
            vm_heads=319.00000000000006,
            thickness_heads=6.510204081632654,
            b_tails=(-0.0002915020000000002-5.0198553394640014e-08j),
            vm_tails=782.0000000000001,
            thickness_tails=15.95918367346939,
            rough_head_tail=3.0,
            rough_preceding_mono=3.0,
            reverse_monolayer=True,
            name='LL rDMPC',
            head_solvent=(6.360408603667384-1.1340361180119794e-07j),
            tail_solvent=(6.360408603667384-1.1340361180119794e-07j),
        )

structure |= LipidLeaflet(
            apm=70.0,
            b_heads=(0.0006007780000000002-2.2512013099632003e-08j),
            vm_heads=319.00000000000006,
            thickness_heads=6.510204081632654,
            b_tails=(-0.0002915020000000002-5.0198553394640014e-08j),
            vm_tails=782.0000000000001,
            thickness_tails=15.95918367346939,
            rough_head_tail=3.0,
            rough_preceding_mono=3.0,
            reverse_monolayer=False,
            name='LL DMPC',
            head_solvent=(-0.5605198489478392-6.185378961109063e-05j),
            tail_solvent=(-0.5605198489478392-6.185378961109063e-05j),
        )
structure |= LipidLeaflet(
            apm=70.0,
            b_heads=(0.0006007780000000002-2.2512013099632003e-08j),
            vm_heads=319.00000000000006,
            thickness_heads=6.510204081632654,
            b_tails=(-0.0002915020000000002-5.0198553394640014e-08j),
            vm_tails=782.0000000000001,
            thickness_tails=15.95918367346939,
            rough_head_tail=3.0,
            rough_preceding_mono=3.0,
            reverse_monolayer=True,
            name='fLL DMPC',
            head_solvent=(-0.5605198489478392-6.185378961109063e-05j),
            tail_solvent=(-0.5605198489478392-6.185378961109063e-05j),
        )

m = SLD((-0.5604469813674761-6.184574861844118e-05j))
structure |= Slab(0.0, m, 3.0, name='H2O')

model = ReflectModel(structure, bkg=0.0)

(The plot shows the orsopy resolved slabs/SLDs as broken lines (equal to the actual LipidLeaflet structure).)

[maxskoda]

Dear @aglavic ,

Many thanks for all the work implementing this. I've not had time to check all the code and behaviour (so I apologise if I didn't quite understand it correctly), but I noticed a slight nomenclature issue specific to the Bilayer class:
the names "outer" and "inner" are slightly confusing if they are intended to represent headgroup and tail materials. Looking at the implementation, outer corresponds to the headgroup material and inner to the tail material, not the two leaflets of the bilayer - is that correct? If so, this will be confusing downstream, since "inner" and "outer" usually refer to the leaflets proximal and distal to the substrate respectively:

Si | SiO2 | Bilayer(DMPC) | solvent 

represents:

 Si | SiO2 | inner head | inner tail | outer head | outer tail | solvent

Also, do I understand correctly that in the consumer application, I can do something like:

blocks = sample.resolve_to_blocks()

for block in blocks:
    if isinstance(block, model_complex.Bilayer):

        head_volume = block.outer.volume
        tail_volume = block.inner.volume
        apm = block.apm

This is great! I will try it at the next possible opportunity.

Thanks again,
Max

[aglavic]

Hi @maxskoda,

Concerning the head / tail nomenclature I might have misunderstood how this is used. I throught it reffered to the material on the inner / outer part of the Bilayer, meaning heads and tails. From what you wrote, I guess you mixed up the order in the representation and it should be:

 Si | SiO2 | inner head | inner tail | outer tail | outer head | solvent

Is that correct?

I would thus suggest to call it heads and tails analogous to the LipidLeaflet.
Should we assume that inner = outer lipid or do it analogous to the hydration with two parameters where the second is None, which defaults to the same material.

Also, is inner/outer a good description? What defines the inner side? The model language can be used for both reflectivity directions, so inner/outer would change place if the model is inverted.

Also, do I understand correctly that in the consumer application, I can do something like:

blocks = sample.resolve_to_blocks()

for block in blocks:
    if isinstance(block, model_complex.Bilayer):

        head_volume = block.outer.volume
        tail_volume = block.inner.volume
        apm = block.apm

Exactly. Have a look at the RefNxResolver class I wrote for refnx examples/resolve_and_plot_model.py to see the details how this could be implemented. It should be a straight-forward mapping to different parameter names with potential unit conversions.

[aglavic]

I've now updated the definition to use heads and tails for the Bilayer. Instead of inner/outer option for different leaflets there is heads_2, tails_2 and hydration_2 as optional second materials defined with respect to the incoming beam.