- Categories
- Cylinder
- polymer_lipid_particle
- polymer_lipid_particle.py
polymer_lipid_particle - polymer_lipid_particle.py
r"""
Definition
----------
This model calculates the scattering from polymer lipid particles (or polymer nanodiscs)s.
Default parameters corresponds to styrene maleic acid (SMA) (2:1) and DMPC.
The model is a circular polydisperse nanodisc model with molecular parametrization through the number of lipids and molecular volumes.
A log normal distribution in the number of lipids was used - directly coded in the model, not through sasview's build-in polydispersity.
b values are the number of electrons, and should correspond to the relevant molecules (default is SMA(2:1) and DMPC), but the b values *should not be fitted*
The ED (or SLD) values are calculated in the model as ED = b/(molecular volume).
The ED is the solvent electron density. Should be adjusted to match buffer (default 1xPBS), but *should not be fitted*.
The model is not on absolute scale, so the form factor is normalized to 1 and the scaling should be set to I(0)
References:
----------
Tang et al (2026) "Structure and stability of styrene-maleic acid lipid particles: the role of polymer-lipid ratio", Journal of Colloid and Interface Science 724: 141054. https://doi.org/10.1016/j.jcis.2026.141054.
Larsen et al (2027) "State-of-the-art analysis of SAXS data of high-density lipoprotein, nanodiscs or polymer-lipid particles – a practical guide" unpublished.
Authorship and Verification
----------------------------
* **Author:** Andreas Haahr Larsen **Date:** 6 November 2026
"""
import numpy as np
from scipy.special import jv
name = "polymer_lipid_particle"
title = "Polymer-lipid particle or polymer nanodisc"
description = """styrene maleic acid nanodiscs 2"""
category = "spheres"
parameters = [["Nlip", "", 130, [20, 500], "", "Number of lipids"],
["rel_poly", "", 0.5, [0.1, 1.0], "", "Relative polydispersity lipids"],
["T", "Å", 8.0, [6, 11], "", "Polymer belt thickness"],
["A_lip", "Å^2", 60.6, [55, 66], "", "Area per lipid (default: DMPC)"],
["sigma_roughness", "Å", 10, [0, 20], "", "Interface roughness"],
["V_tail", "Å^3", 743, [670, 820], "", "Volume of lipid tail (default: DM of DMPC)"],
["V_head", "Å^3", 327, [300, 360], "", "Volume of lipid head (default: PC of DMPC)"],
["V_polymer", "Å^3", 461, [360, 560], "", "Volume of polymer unit (default: SMA2000 (2:1))"],
["b_tail", "#electrons", 210, [0, np.inf], "", "Electrons per lipid tail (default: DM of DMPC)"],
["b_head", "#electrons", 164, [0, np.inf], "", "Electrons per lipid head (default: PC of DMPC)"],
["b_polymer", "#electrons", 171, [0, np.inf], "", "Electrons per polymer subunit (default: SMA2000 (2:1))"],
["ED_solvent", "#electrons/Å^3", 0.336, [0.334,0.338], "", "Electron density solvent (default: 1xPBS buffer)"]
]
def Iq(q,Nlip,rel_poly,T,A_lip,sigma_roughness,V_tail,V_head,V_polymer,b_tail,b_head,b_polymer,ED_solvent):
"""
Description
"""
# electron densities (scattering length densities obtained by multiplying with classic electron radius 2.82 fm)
V_lip = V_tail+V_head
SLD_H2O = ED_solvent # 0.336 e/Å^3 for the buffer instead of 0.334 for pure H2O
SLD_tail = b_tail/V_tail
SLD_head = b_head/V_head
SLD_polymer = b_polymer/V_polymer
# dSLDs
dSLD_head = SLD_head - SLD_H2O
dSLD_tail = SLD_tail - SLD_H2O
dSLD_polymer = SLD_polymer - SLD_H2O
# renaming of dSLDs in terms of rho (written as p here in the code)
p_core = dSLD_tail
p_shell = dSLD_head
p_rim = dSLD_polymer
# bilayer height
H_core = 2*V_tail/A_lip # height of lipid core, bilayer, therefore factor 2
H_lip = 2*V_lip/A_lip # total lipid height
H_rim = H_core # hydrophobic matching
# polydispersity in number of lipids (and effectively in number of SMAs)
distribution = "lognormal" # options: "normal" or "lognormal"
N_lip_mu = Nlip
# Compute distribution parameters (different polydispersity options)
if distribution in ["normal","Normal","Gauss","gauss","Gaussian","gaussian"]:
N_lip_sigma = rel_poly * N_lip_mu
N_lip_min = max([N_lip_mu - 3 * N_lip_sigma, 10]) # min 10 Nlips
N_lip_max = N_lip_mu + (N_lip_mu - N_lip_min)
weights = lambda N: np.exp(-0.5 * ((N_lip_mu - N) / N_lip_sigma) ** 2)
steps = 20
elif distribution in ["lognormal","log-normal","LogNormal","Log-Normal","loggauss"]:
sigma_ln = np.sqrt(np.log(1 + rel_poly**2))
mu_ln = np.log(N_lip_mu) - 0.5 * sigma_ln**2
N_lip_min = np.exp(mu_ln - 2.33 * sigma_ln) # 1% quantile
N_lip_max = np.exp(mu_ln + 2.33 * sigma_ln) # 99% quantile
weights = lambda N: lognormal_pdf(N, mu_ln, sigma_ln)
steps = 20
else:
raise ValueError("Unsupported size distribution type")
if rel_poly:
N_lip_array = np.linspace(N_lip_min, N_lip_max, steps)
weights_array = np.array([weights(N) for N in N_lip_array])
weights_norm = weights_array.sum()
weights_array /= weights_norm
N_lip_theory = np.linspace(0.0001,int(N_lip_max+100),int(N_lip_max*10))
weights_theory = np.array([weights(N) for N in N_lip_theory])
weights_theory /= weights_norm
else:
N_lip_array,weights_array = np.array([N_lip_mu]),[1]
N_lip_theory,weights_theory = np.array([N_lip_mu*0.9,N_lip_mu*0.999,N_lip_mu,N_lip_mu*1.001,N_lip_mu*1.1]),[0,0,1,0,0]
A_core = N_lip_array/2.0 * A_lip
R_core = np.sqrt(A_core/np.pi)
A = P_nanodisc_vector(q, R_core, H_core, H_lip, H_rim, T, p_core, p_shell, p_rim, weights_poly=weights_array)
A /= A[0] # normalization
I_total = A * np.exp(-0.5*(q*sigma_roughness)**2)
return I_total
Iq.vectorized = True # Iq accepts an array of q values
##################################################################################
### support functions
##################################################################################
def lognormal_pdf(x, mu, sigma):
# return exp(-((log(x) - mu) ** 2) / (2 * sigma ** 2)) / (x * sigma * sqrt(2 * pi))
f = np.exp(-((np.log(x) - mu) ** 2) / (2 * sigma ** 2)) / x
norm = sigma * np.sqrt(2 * np.pi)
return f/norm
def P_nanodisc_vector(q, R_core, H_core, H_lip, H_rim, T_rim, p_core, p_shell, p_rim, weights_poly=None):
"""
Vectorized form factor: nanodisc
- Integrates over orientation (alpha) using Gaussian quadrature
- Optionally integrates over polydispersity (R_core_array)
Parameters
----------
q : array: Scattering vector magnitudes
R_core : Array: of core radii (for polydispersity). Can be length 1 for no polydispersity
H_core, H_lip, H_rim, T_rim : float: Geometric parameters
p_core, p_shell, p_rim : float: Excess scattering length densities
N_alpha : int: Number of Gaussian quadrature points for orientation integration
weights_poly : array: weights for polydispersity. If None, no polydispersity weighting applied
Returns
-------
A : array: Scattering amplitude
"""
## sample alpha and calculate their weights
alpha, weights_alpha = gauss_legendre_alpha(N_alpha=12, alpha_min=0, alpha_max=np.pi/2)
## initialize arrays for vectorized calculations
q_b = q[:, None, None] # (len_q,1,1)
alpha_b = alpha[None, :, None] # (1,N_alpha,1)
weights_alpha_b = weights_alpha[None, :, None] # (1, N_alpha, 1)
R_b = R_core[None, None, :] # (1,1,len_R_core)
weights_poly = np.asarray(weights_poly)[None, :] # (1, len_R_core)
## compute psi1 for all q, alpha, R_core
A_nd = psi_nanodisc(q_b, R_b, H_core, H_lip, H_rim, T_rim, p_core, p_shell, p_rim, alpha_b)[1] #(len_q, ??,??)
## integrate over alpha
A_alpha = np.sum(np.abs(A_nd)**2 * weights_alpha_b, axis=1) # (len_q, len_weights_poly), integrate out orientation dimension (alpha)
## integrate over polydispersity
A = np.sum(A_alpha * weights_poly, axis=1) # (len_q), integrate out polydispersity dimension
return A
def P_nanodisc(q,R_core,H_core,H_lip,H_rim,T_rim,p_core,p_shell,p_rim):
"""
form factor: nanodisc
R_core: radius of core
H_core: height of core
H_lip: total height of bilayer (core+shell)
H_rim: height of rim/belt
T_rim: thickness of rim/belt
p_core: excess scattering length density (DeltaSLD) for the core, i.e. lipid tail groups
p_shell: DeltaSLD for shell, i.e. lipid head groups
p_rim: DeltaSLD for rim/belt
"""
"""
alpha_min,alpha_max,N_alpha = 0,pi/2,7
alpha = np.linspace(alpha_min,alpha_max,N_alpha)
d_alpha = alpha_max/(N_alpha-1)
sin_alpha = sin(alpha)
P = np.zeros(len(q))
A = np.zeros(len(q))
for (a,sin_a) in zip(alpha,sin_alpha):
psi0,psi1 = psi_nanodisc(q,R_core,H_core,H_lip,H_rim,T_rim,p_core,p_shell,p_rim,a)
P += abs(psi0)**2 * sin_a
A += abs(psi1)**2 * sin_a
P *= d_alpha
A *= d_alpha
"""
N_alpha = 7
x, w = np.polynomial.legendre.leggauss(N_alpha)
alpha = 0.25 * np.pi * (x + 1) # map [-1,1] → [0,pi/2]
weights = 0.25 * np.pi * w * np.sin(alpha) # include Jacobian + sin(alpha)
psi0,psi1 = psi_nanodisc(q[:, None], R_core, H_core, H_lip, H_rim, T_rim, p_core, p_shell, p_rim, alpha[None, :])
P = np.sum(np.abs(psi0)**2 * weights, axis=1)
A = np.sum(np.abs(psi1)**2 * weights, axis=1)
return P,A
def gauss_legendre_alpha(N_alpha=15, alpha_min=0, alpha_max=np.pi/2):
"""
Generate Gauss–Legendre quadrature nodes and weights
for orientation averaging over [alpha_min, alpha_max].
Parameters
----------
N_alpha : int
Number of quadrature points (default 15)
alpha_min : float
Lower bound of alpha integration (default 0)
alpha_max : float
Upper bound of alpha integration (default pi/2)
Returns
-------
alpha : ndarray
Quadrature nodes mapped to [alpha_min, alpha_max]
weights_alpha : ndarray
Quadrature weights including Jacobian sin(alpha)
"""
# Gauss–Legendre nodes and weights on [-1,1]
x, w = np.polynomial.legendre.leggauss(N_alpha)
# Map to [alpha_min, alpha_max]
alpha = 0.5 * (alpha_max - alpha_min) * (x + 1) + alpha_min
# Scale weights and include Jacobian sin(alpha)
weights_alpha = 0.5 * (alpha_max - alpha_min) * w * np.sin(alpha)
return alpha, weights_alpha
def psi_nanodisc(q,R_core,H_core,H_lip,H_rim,T_rim,p_core,p_shell,p_rim,a):
"""
form factor amplitude: nanodisc
R_core: radius core
H_core: height core (2xlipid tail group)
H_lip: total height (core and shells)
H_rim: height of rim/belt (protein or polymer)
T_rim: thickness of rim/belt (protein or polymer)
p_core: excess scattering length density (DeltaSLD) for the core, i.e. lipid tail groups
p_shell: DeltaSLD for shell, i.e. lipid head groups
p_rim: DeltaSLD for rim/belt
a: alpha (orientation parameter) - can be array
"""
# lipid bilayer core-shell center
V_core = V_cyl(R_core,H_core)
psi_core = psi_cyl(q,R_core,H_core,a)
V_lip = V_cyl(R_core,H_lip)
psi_lip = psi_cyl(q,R_core,H_lip,a)
A_core = p_core * V_core * psi_core
norm_core = p_core * V_core
A_shell = p_shell * ( V_lip * psi_lip - V_core * psi_core )
norm_shell = p_shell * ( V_lip - V_core )
# belt or rim of polymers/proteins/peptides
V_outer = V_cyl(R_core+T_rim,H_rim)
psi_outer = psi_cyl(q,R_core+T_rim,H_rim,a)
V_inner = V_cyl(R_core,H_rim)
psi_inner = psi_cyl(q,R_core,H_rim,a)
A_rim = p_rim * ( V_outer * psi_outer - V_inner * psi_inner )
norm_rim = p_rim * ( V_outer - V_inner )
# collect all terms
A_nanodisc = A_core + A_shell + A_rim
norm_nanodisc = norm_core + norm_shell + norm_rim
psi_nanodisc = A_nanodisc/norm_nanodisc
return psi_nanodisc,A_nanodisc
def V_cyl(R,L):
"""
Volume of a cylinder with radius R and length L
Supports scalar or vector R. q and a can also be vectors.
"""
return np.pi*R**2*L
def psi_cyl(q,R,L,a):
"""
form factor amplitude: cylinder with radius R and length L
depends on view angle a (for alpha)
reference: Pedersen1997
"""
x = q*R*np.sin(a)
z = q*L*np.cos(a)/2.0
psi=2.0*bessj1c(x)*sinc(z)
return psi
def sinc(x):
"""
function for calculating sinc = sin(x)/x
numpy.sinc is defined as sinc(x) = sin(pi*x)/(pi*x)
"""
return np.sinc(x/np.pi)
def bessj1(x):
""""
Bessel function of the first kind of first order
"""
return jv(1,x)
def bessj1c(x):
"""Calculate J1(x)/x for N-dimensional arrays"""
x = np.asarray(x)
y = np.full_like(x, 0.5, dtype=float)
mask = x != 0
y[mask] = bessj1(x[mask]) / x[mask]
return y
##################################################################################
### TEST:
##################################################################################
if __name__ == '__main__':
import time
print(f'name: {name}')
print(f'title: {title}')
print(f'description: {description}')
for p in parameters: assert len(p[3]) == 2, f"Typo in {p[0]} bounds: expected [min, max], got {p[3]}"
start_time = time.time()
q = np.linspace(0.0001,1.2,5000)
p0 = []
for p in parameters:
p_name,p_units,p_guess,p_bounds,p_empty,p_description = p
p0.append(p_guess)
I = Iq(q,*p0) + 1e-3
end_time = time.time()-start_time
print(f'Finished, time: {end_time:.2f} s')
if 1:
import matplotlib.pyplot as plt
plt.figure()
plt.plot(q,I,label='total intensity')
plt.xscale('log')
plt.yscale('log')
plt.legend ()
plt.show()
Back to Model
Download