Since the solution of a differential equation on the lattice is about ten thousand times faster than a single simulation we compare the results obtained with the solution of the elastic model to be able to explore the large parameter space.
We refer to the continuum model for bilayer membranes defined by [Westet al.(2009);
Branningan and Brown (2006, 2007)] where the membrane is described as two thick coupled elastic sheets. Each elastic sheet has an energy cost for bending and stretching and the elastic coupling is described as an harmonic potential. The free energy of the
model is: kbenis the bending rigidity,γ0
controls the asymmetry between the bilayer (it is connected with the spontaneous curvature,kel
the elastic coupling.
continuum model CHAPTER 6. INCLUSION INTERACTIONS The term γ0 is connected with the spontaneous curvature, c0, with this relationship [Branningan and Brown (2006)]:
γ0 := 4kben
c0−∂c0
∂AA
(6.4) The Hamiltonian (eq:6.3) recalls the Hamiltonian already introduced in (eq:4.32) where we have excluded the protrusion contributions (dependent on the coupling termsγλ, kλ) that can not be described in continuum models. In (4.32) we have definedh(x, y) as the height of the membrane at the point (x, y) and d(x, y) as the thickness of the membrane at the same point. We have changed the coordinate system and disregarded the fluctuations over the averageh(x, y). In (6.3) we considerds(x, y) as the position of the intersection between the hydrophobic and the hydrophilic beads of the upper leaflet (the lower leaflets is symmetric to the upper one). The transformation between d(x, y) andds(x, y) is simply:
ds(x, y) =d(x, y)−d¯ (6.5)
This simple change allows us to describe the membrane modification by a single func-tion, ds(x, y), and ds(x, y) = 0 corresponds to the ground state of the energy (6.3).
Every modification of the surface ds(x, y) from the ground state gives a energetic contribution in terms of bending, stretching and elastic coupling.
Figure 6.2: lhs) Sketch of the numerical solution of the Euler-Lagrange equations in 2d. The bonds represent the connections between the lattice points. The blue and red points fix the boundary conditions. The violet points represent the thickness profile calculated from the simulation on the same grid discretisation as the numerical solution.
Euler-Lagrange equations
We calculate theEuler-Lagrangeequation from the free energy expressed in (see eq:6.3) and obtain [Westet al.(2009)]:
δF
δds = 0 (kben∇4+γ0∇2+kel)ds(x, y) = 0 (6.6) We solve the differential equation (eq:6.7) in different steps. We consider first the Cartesian one dimensional problem.
Dxds(x) = 0 Dx:=kben∇4x+γ0∇2x+kel ∇2x:= d2
dx2 (6.7)
To calculate the solution of the linear differential equation we search for a solution of the kind exp(λx) and write the characteristic equation:
λ±± are the roots of the characteristic equation
kbenλ4+γ0λ2+kel = 0 λ±±=± s
−γ0±p
γ02−2kbenkel
2kben (6.8)
which gives the solution:
ds(x) =c1eλ++x+c2e−λ++x+c3eλ+−x+c4e−λ+−x (6.9)
continuum model CHAPTER 6. INCLUSION INTERACTIONS We see that the following boundary conditions:
ds(rp) =hp d0s(rp) = tanθc ds(Lx) = 0 d0s(Lx) = 0 (6.10) are sufficient to determine the parametersci.
Second we calculate the radial solution where the Laplacian operator turns into:
Drds(r) = 0 Dr:=kben∇4r+γ0∇2r+kel ∇2r:= 1
r∂r(r∂r) (6.11)
The radial solution of the equation (6.7) is [Westet al.(2009)]: ∇2r=1r∂r(r∂r) in cylindrical coordinates,J0(r) andY0(r) are the 0th-orderBesselfunctions of the first and second kind.
ds(r) =c1J0(λ++r) +c2Y0(λ++r) +c3J0(λ+−r) +c4Y0(λ+−r) (6.12) where the λ coefficients are the same calculated for the one dimensional Cartesian problem (eq:6.9) and we use the same boundary conditions (eq:6.10).
Third we define the differential operator in terms of partial derivatives in Cartesian coordinates:
Dcds(x, y) = 0 Dc:=kben∇4c+γ0∇2c+kel ∇2c :=
∂x2 0 0 ∂y2
(6.13) The solution of partial differential equations is a complex problem that uses theLax -Milgramtheorem to show the existence and uniqueness of a weak solution to a given boundary problem [Babuˇska (1971)]. To avoid the complexity of treating the mathe-matical problem we test the solution of the discretized two dimensional problem with the analytical radial solution. We create a two dimensional square lattice, we dis-cretized the operators and we put the following boundary conditions:
{ds(x, y) =hp|x2+y2=rp} {d0s(x, y) = tanθc|x2+y2=rp}
{ds(x, y) = 0|x=Lx∨y=Ly} {d0s(x, y) = 0|x=Lx∨y=Ly} (6.14) We calculate the numerical solution of the Euler-Lagrange equations for a square lattice in two dimensions (see fig:6.2) using finite differences (see appendix).
Figure 6.3: lhs) thickness profile in Cartesian
coordinates. rhs) comparison between the thickness profiles.
The green and the purple lines represent the radial average from the Cartesian thickness profiles form the simulations (green line) and from the 2d numerical solution of (eq:6.15) (red line). The red line is obtained by fitting with the simulation data with the functionAe−x/λsin(f x)/x, the blue line shows the numerical solution in cylindrical coordinates.
To map the parameter of the differential equation into the simulated system we rescale the constants:
∇4+ γ0 kben
∇2+ kel kben
ds(x, y) = 0 (6.15)
simulation description CHAPTER 6. INCLUSION INTERACTIONS We refer to a lipid bilayer using the values calculated in the previous chapter (see table:4.3) and use the following values:
∇4+ 0.5∇2+ 0.57
ds(x, y) = 0 (6.16)
If the peptide sits in the middle of the lattice we recover radial symmetry around the peptide center. We reduce the two dimensional solution into the radial one integrating over the angle and we superimpose the two solutions on the same plot (see fig:6.3).
The thickness profile obtained from the simulations is the result of the radial average of the density profile. We have rescaled the density profile by the ratio ¯d/ρ¯and shifted the value of the bulk thickness, ¯d, to the ground line:
d(r) = ¯d ρA(r)
¯ ρ −1
(6.17) The thickness profile obtained from the numerical calculation is the radial average profile of the height of the grid points around the inclusion.
The agreement between the analytical solution, the reduced two dimensional solu-tion and the thickness profile obtained from the simulasolu-tions (see fig:6.3) furnishes a satisfying argument about the reliability of the two dimensional solution.
We want to express the solution using a physical analogy with the damped harmonic oscillator in far field approximation in cylindrical coordinates:
Aamplitude,ξcorrelation length andf the oscillation wave vector.
ds(r) =Ae−r/ξsin(f r)
r (6.18)
Fitting the solution of the numerical calculation or the simulation we can extract the correlation length and the oscillation wave vector of the protein lipid interaction:
ξ= 2.25[∆L] andf = 0.88[∆L−1] (see fig:6.3).