Kinetics of self-assembly of inclusions due to lipid membrane thickness interactions

Self-assembly of proteins on lipid membranes underlies many important processes in cell biology, such as, exo- and endo-cytosis, assembly of viruses, etc. An attractive force that can cause self-assembly is mediated by membrane thickness interactions between proteins. The free energy profile associated with this attractive force is a result of the overlap of thickness deformation fields around the proteins. The thickness deformation field around proteins of various shapes can be calculated from the solution of a boundary value problem and is relatively well understood. Yet, the time scales over which self-assembly occurs has not been explored. In this paper we compute this time scale as a function of the initial distance between two inclusions by viewing their coalescence as a first passage time problem. The first passage time is computed using both Langevin dynamics and a partial differential equation, and both methods are found to be in excellent agreement. Inclusions of three different shapes are studied and it is found that for two inclusions separated by about hundred nanometers the time to coalescence is hundreds of milliseconds irrespective of shape. Our Langevin dynamics simulation of self-assembly required an efficient computation of the interaction energy of inclusions which was accomplished using a finite difference technique. The interaction energy profiles obtained using this numerical technique were in excellent agreement with those from a previously proposed semi-analytical method based on Fourier-Bessel series. The computational strategies described in this paper could potentially lead to efficient methods to explore the kinetics of self-assembly of proteins on lipid membranes. Author summary Self-assembly of proteins on lipid membranes occurs during exo- and endo-cytosis and also when viruses exit an infected cell. The forces mediating self-assembly of inclusions on membranes have therefore been of long standing interest. However, the kinetics of self-assembly has received much less attention. As a first step in discerning the kinetics, we examine the time to coalescence of two inclusions on a membrane as a function of the distance separating them. We use both Langevin dynamics simulations and a partial differential equation to compute this time scale. We predict that the time to coalescence is on the scale of hundreds of milliseconds for two inclusions separated by about hundred nanometers. The deformation moduli of the lipid membrane and the membrane tension can affect this time scale.

This paper is arranged as follows. First, we quantify the interaction energy profile of 54 hexagonal, rod-and star-shaped inclusions 1 . We show that our finite difference 55 numerical method for computing energies agrees very well with analytical formulae 56 (using Fourier-Bessel series) in most cases. After computing the interaction energies, we 57 solve first-passage time problems to find the time scales over which two inclusions 58 coalesce due to attractive interactions. We use both Langevin dynamics and the 59 Fokker-Planck equation to solve first passage time problems and study both isotropic 60 and anisotropic problems with reflecting/absorbing boundary conditions. Finally, we 61 summarize our results in the discussion and conclusion sections and point to various 62 enrichments that can be implemented following our earlier work [26].  [11] R 1 (θ 1 ) shape function for the centered inclusion nm R 2 (θ 2 ) shape function for the moving inclusion nm θ the angle between two inclusions and horizontal line (see Fig 2(b)) radian/degree u thickness deformation nm u i thickness deformation at node i nm A e the area of the triangular element at node i nm 2 R 1 radius of the inner boundary w.r.t (r 1 , θ 1 ) nm R 2 radius of the outer boundary w.r.t (r 2 , θ 2 ) nm u b the vector of all nodes determined by Eq (14)  We consider a circular lipid membrane with radius R 2 and two inclusions embedded in 66 it. Our first goal is to compute the energy landscape seen by an inclusion interacting 67 with another inclusion on a flat membrane. The interactions between the inclusions are 68 a result of the overlap of membrane thickness deformation fields in their vicinity. The 69 interaction energy will be computed by considering two inclusions, one fixed and the 70 other moving as shown in Fig 2(a). The coordinate frame at the fixed inclusion (blue) 71 denoted as inclusion 1 (r 1 , θ 1 ) is set to be the default one. Assume that the moving 72 inclusion (purple) denoted as inclusion 2 initially stays in the same orientation as 73 inclusion 1 (see Fig 2(a)). To keep the analysis simple, when an inclusion moves we do 74 not consider its rotational diffusion. As inclusion 2 moves from the initial configuration 75 to the green spot and forms an angle θ with the horizontal line (see Fig 2(b)), the 76 September 17, 2020 4/23 energy of the system can be computed by rotating both inclusions anticlockwise by 77 angle θ from the initial configuration (see Fig 2(c)). This interaction energy will enter 78 our analysis of the kinetics of the moving inclusion due to Brownian motion.
In order to efficiently apply the boundary conditions, we rewriteū 2 as a function of 92 r 1 , θ 1 , r andū 1 as a function of r 2 , θ 2 , r, We consider the 94 following type of boundary conditions (see Fig 1( We can solve for the 4(2N + 1) coefficients 96 A ± 1,0 , A ± 2,0 , A ± 1,n , A ± 2,n , B ± 1,n , B ± 2,n , n = 1, 2, · · · , N because Eq (9)-(10) result in a linear 97 system. This determines the full deformation field due to the overlap of the 98 deformations caused by both inclusions. The next step is to compute the energy φ(r) 99 due to this deformation field. Note that the angular dependence of φ(r) appears in the 100 shape functions of two inclusions, R 1 , R 2 .

101
Using the divergence theorem, the total energy expression in Eq (1) can be converted 102 to the sum of line integrals over the boundary terms, i.e. φ = φ 1 + φ 2 with φ i given by, 103 From the first line to the second line we assume the line integral along the outer 104 boundary is a constant (which works out to 0 as R 2 → ∞ and F → 0) w.r.t r and put it 105 into the G 1 term (both G 0 and G 1 are constants). The energy φ(r) can be computed 106 relatively efficiently using this technique. This is important since φ(r) must be 107 computed repeatedly as inclusion 2 moves and r changes due to Brownian motion when 108 we solve the first passage time problem. We will also need the forces acting on inclusion 109 2 in our analysis later. Eq (11) gives an expression to compute the force analytically, 110 which in the special case of an isotropic φ(r) (i.e., no angular dependence) works out to 111 When there is only one circular inclusion in the membrane, the thickness deformation 112 field in Eq (5) has a closed form solution [7] which can be compared to the simulation 113 result of Klingelhoefer et al. [28] who studied radial bilayer thickness profiles for the Gα 114 nanopore (among many others). We used the same parameters and boundary conditions 115 as they did: a = 34.19Å, U 1 = 0.81Å, U 1 = 0.7, R 1 = 10Å and fit their curves by  Finite difference method based on refined grid 119 The above analysis gives us a semi-analytical technique to compute φ(r). However, not 120 all problems can be solved analytically, so a computational method is needed to 121 estimate interaction energies. Fortunately, Eq (1) can be minimized using a finite 122 difference method. We discretize the membrane using equilateral triangle elements as 123 in [29,30]. The thickness deformation at node i is denoted by u i . Since the thickness 124 deformation u changes rapidly in the neighborhood of the two inclusions, we used 125 refined grids in a region containing the two inclusions and non-refined grids in the  Using methods similar to those in [29,30] the energy is first written in a discrete 129 form and then the thickness deformation field that minimizes this energy is computed. 130 Finally, the minimizer is plugged back into the energy expression, and we have, where 1 u is a column vector of size len(u) with all entries 1 and λ = 2F Ae a . It can be 132 shown that the boundary condition in [22] can be written in the discrete form, k pairs along the boundary (see Fig 4(a)). (15) Note that u i , u j are given in Eq (14) and thus u k can be sloved from Eq (15) 134 immediately. We also assume the inclusions are stiff such that for all nodes i inside the 135 inclusions u i are a constant (equal to the value of those at boundary nodes). Hence,

136
Eq (13) can be rewritten as, Taking (16) is 138 minimized where 1 a is a column vector of size len(u a ) with all entries 1. Then, we can 139 write the minimized total energy as, Applications to hexagon, rod and star shaped inclusions 141 We now focus on the interaction of two hexagon shaped inclusions on a lipid membrane 142 which has a rotational periodicity of π/3. We use Eq (17) derived from our numerical 143 method and Eq (11) derived from the analytical method, to compute the interaction 144 energy of two inclusions with distance r and then make comparisons. Agreement of the 145 results from the two methods is expected. As shown in Fig 5(a), the energy computed 146 using the analytical method for two inclusions separated by distance r in two different 147 orientations (shown in the inset differing by a rotation of π/6) are almost the same.

148
Hence, we can simplify our model and consider the energy landscape generated by two 149 hexagon inclusions as being almost isotropic (insensitive to rotation). In Fig 5(b), we fix 150 the shape of two hexagons (see inset of the figure), and show that as the number of grid 151 elements per side n increases, the match between the energy computed from the 152 numerical method and analytical method gets better, justifying our numerical approach 153 of using refined grids near the inclusions. From Fig 5(c) we learn that as applied tension 154 increases, the attraction at small separations (around r = R 1 = 7 nm) becomes weaker 155 but the repulsive force at around 9 − 10nm also becomes weaker. Please check those to 156 be sure.    shows that as applied tension increases, the force becomes weaker at short separations, 167 which implies that elastic interactions could be weakened by strong applied tension.

168
Physically, this is reasonable, since high tension will tend to make the membrane flatter 169 so that the thickness is more uniform everywhere. Our main goal in this paper is to study the kinetics of an inclusion diffusing in an 186 energy landscape resulting from elastic interactions with another inclusion. Efficient 187 methods to compute the energy landscape developed above are a pre-requisite to this 188 exercise. We will now use these methods to solve first passage time problems. 189 We consider a circular membrane of size R 2 = 125 nm with a circular inclusion of 190 size 2.5 nm fixed at the center. Another circular inclusion of the same size is diffusing 191 around driven by stochastic forces. Recall from the energy landscape that there are 192 attractive interactions between inclusions when the separations are small. Hence, if the 193 moving inclusion comes close enough to the static one at the center then it will be 194 strongly attracted. Therefore, we assume that at R 1 = 7 nm there is an absorbing wall 195 at which the moving inclusion will disappear by being attracted towards the center. We 196 assume that at R 2 = 125 nm (far away) there is a reflecting wall where the moving 197 inclusion will be bounced back. Note that problems in which both boundaries are 198 absorbing were solved elsewhere [26]. The exercise we will perform now is as follows. 199 We place the second inclusion randomly on a circle of radius r = y at time t = 0 and let 200 it diffuse around. At some time t = T in when the inclusion hits the inner boundary for 201 the first time we stop it from diffusing and record T in . We repeat this experiment a 202 large number of times and record T in for each repetition. The mean value of T in is the 203 mean first passage time T 1 . Our goal is to find T 1 (y) as a function of the initial 204 condition r = y. This can be done analytically or through a Langevin dynamics 205 simulation. We will use both methods in the following.

206
To estimate T 1 (y) analytically we first need to compute survival probabilities. Let p 207 be the probability density (for finding the inclusion) at position r and angle θ given 208 initial condition r = y, θ = α and P (r, t|y) = with Dirichlet boundary condition at the inner boundary and Robin boundary condition 213 at the outer boundary [31], In the above D is a diffusion coefficient of the inclusion in the lipid membrane and ν is 215 a drag coefficient which are related by the Nernst-Einstein relation νD = k B T [26]. Let 216 S(y, t) be the survival probability, Then, we can get the first passage time density, The existence of the first moment of P (r, t|y) with respect to time t can be shown from 219 the fundamental solution constructed by Itô in [32]. Then, tP (r, t|y) → 0 as t → ∞.

220
Accordingly, the first passage time T 1 (y) can be derived from Eq (21), where g 1 is defined by, Theorem 1: The ODE for T 1 (y) with a reflecting wall at the outer boundary and an 223 absorbing wall at the inner boundary is with boundary conditions, Proof: See Proof of Theorem 1.

226
Next, we describe how to estimate T 1 (y) using Langevin dynamics simulations. The 227 overdamped version of the Langevin equation in an isotropic setting is given by [26], where i represents two perpendicular directions of the motion. ν is the translational 229 drag coefficient of a circular inclusion given by the Saffman-Delbrck model [33]. represents the stochastic force along direction i. We initially put the moving particle 232 somewhere at r = y, and choose a time step dt that ensures convergence of the Lagenvin 233 dynamics simulation. Then, for each time step dt, we perform the calculation in 234 Eq (26), updating the position of the moving inclusion. We record the time at which 235 the moving particle hits the absorbing wall at R 1 . We run 8000 simulations and then 236 take an average to estimate the first passage time. For more details the readers are 237 referred to [26]. Fig 5(a) and Fig 8(a) show that the φ(r) for hexagon and star 238 inclusions can be regarded as nearly isotropic. We use Eq (24) to numerically solve for 239 the first passage time and compare the results obtained from the two methods.  For two non-circular inclusions, the corresponding Fokker-Planck equation is a partial 254 differential equation of parabolic type [26], Accordingly, we need to redefine the first passage time in Eq (22) which is now given by, 256 Theorem 2: The PDE for T 1 (y, α) with a reflecting wall at the outer boundary and an 259 absorbing wall at the inner boundary is given below, Proof: See Proof of Theorem 2.

262
The overdamped Langevin equation in an anisotropic setting is given by [26], where i represents two perpendicular directions of the motion and ν ij and ξ i are the 264 diffusion tensor and random force tensor (for more details see [26]). Fig 6(a) shows that 265 the interaction energy φ(r) for rod shaped inclusions depends on θ (it is anisotropic). In 266 the Langevin dynamics calculations, for each initial position y we use Eq (33) to run The good agreement between the first passage time solved from the PDE in Eq (31) 272 and estimated by Langevin equation once again shows that our methods are work well. 273 As shown in Fig 11, as the initial angle increases from 0 • to 90 • , the first passage time 274 decreases at small separations but increases at large separations. This can be similarly 275 explained by the fact that stronger attractive force near R 1 pulls the moving particle to 276 be absorbed faster from smaller initial separations while stronger repulsive force around 277 12 − 16 nm leads to a larger first passage time when the particle is initially located at a 278 large distance.    This paper has two major parts. In the first part we use a finite difference method to 286 compute the interaction energy of two inclusions due to membrane thickness 287 deformations. In the second part we use the computed energy landscape to solve first 288 passage time problems. Our method to compute energies is different from the analytical 289 method in [7,8] which uses perturbation theory to study thickness mediated interactions 290 between two anisotropic inclusions; we implement an approach to compute the energy 291 using the divergence theorem which is more general and can deal with strongly 292 anisotropic inclusions. The advantage of analytical methods in both [7,8] and this work 293 is that they can compute the energy accurately at small applied tension F if enough 294 terms in the Fourier-Bessel series are used. However, it is time consuming to compute 295 the coefficients in the Fourier-Bessel series and this becomes computationally infeasible 296 when the inclusions are strongly anisotropic. On the other hand, our numerical method 297 is able to handle arbitrary values of F and can efficiently compute the interaction 298 energy of two inclusions for different separations r given a fixed set of parameters 299 (K b , K t , a etc.) which are stored in a pre-calculated stiffness matrix.

300
In the second part of the paper we compute the time to coalescence of two inclusions 301 of various shapes as a function of the distance separating them. We use both Langevin 302 dynamics and a PDE to arrive at our estimates. For two inclusions separated by about 303 125 nm we predict that the time to coalescence is hundreds of milliseconds irrespective 304 of the shape of the inclusion. The time to coalescence with only membrane bending 305 interactions was of similar magnitude as shown in [26]. The order of magnitude of the 306 time to coalescence is the same even though the attractive force due to membrane 307 thickness interactions is stronger than that due to membrane bending interactions in [26] 308 at small separations. The reason is that even with membrane thickness interactions the 309 attractive force decays to zero quickly and Brownian motion dominates the kinetics of 310 the moving particle in most regions, just as in [26]. Therefore, at small separations the 311 first passage time with thickness mediated interactions is smaller than that with 312 out-of-plane bending interactions, but is not very different at large separations.

314
In this paper we have analyzed the temporal self-assembly of inclusions due to 315 interactions mediated by membrane thickness variations. A major accomplishment of 316 the paper is to show that the results from Langevin dynamics simulations agree well 317 with those obtained from a PDE for the first passage time. The approach based on the 318 PDE is much faster than the Langevin dynamics simulation and could open new ways 319 to study the process of self-assembly. This is a step beyond earlier studies which focused 320 on the energy landscape of clusters of proteins, but did not look into kinetics. Some 321 papers based on molecular simulation did consider the temporal process, but to the best 322 of our knowledge most did not reach the time scales calculated in this paper. We close 323 this paper by mentioning some effects that we did not consider. First, hydrodynamic 324 interactions between inclusions (based on the Oseen tensor) were shown to speed up 325 self-assembly in [26] and they are expected to have a similar effect here. Second, the 326 temporal behavior of a cluster of inclusions are not studied in this paper due to 327 limitations of computational power, but we expect the overall behavior to be similar to 328 the clusters studied in our earlier work [26]. Third, only a limited set of inclusion shapes 329 are considered in this paper, but it is found that the time to coalescence does not 330 depend strongly on shape. We leave it to future work to add these refinements and 331 extend this type of analysis to important functional proteins such as ion-channels [22]. 332

333
Proof of Theorem 1 Following techniques in [26,31,34], we integrate Eq (18) for P 334 over all t ≥ 0, where 1 r δ(r − y) is the initial condition and the second order linear differential operator 336 with domain Using the method in [26], we can get the adjoint operator L * r which satisfies and the inner product is defined as, Proofs for the existence of the solutions of second order inhomogeneous linear ordinary 343 differential equation are well known. Hence, we can find a u 0 ∈ C 2 ([R 1 , R 2 ]) s.t.
Using Eq (38), we can derive Eq (24), a second order ODE for T 1 (y). The boundary 346 condition of T 1 (y) at the absorbing wall is straightforward [31,34]: T 1 (R 1 ) = 0. For the 347 boundary condition at the reflecting wall, we appeal to the Langevin equation in 348 Eq (26). If the particle sits at position R 2 , decomposing the overdamped Langevin 349 equation [26] into radial direction and angular direction, we have, After time dt, the particle can only move to R 2 + dy(dy < 0) along the radial direction 351 because of the reflecting wall at R 2 . The motion along the angular direction can be 352 neglected because T 1 (y) does not have dependence on angular direction. Note that dy is 353 a random variable depending on ξ y and dt with constraint where we used mean value theorem twice to reach to Eq (45) with 356 R 2 + dy < R 2 + η dy < R 2 + β dy < R 2 . Note that β dy depends on η dy and thus depends 357 on dy. C 2 (dt) is the value to satisfy R 2 + dy = R 1 for given dt and ξ y . C 1 (dt) is the 358 scaling factor such that the integral of probability density equals 1: After some re-arrangements and 360 dividing by dt on both sides, As dt → 0, C 1 → 2, C 2 → −∞. Note that η dy dy < 1 and we have |T 1 (R 2 + β dy )| < M 362 for some M because T 1 is C 2 . Then if we take t → ∞ the third term in RHS of Eq (46) 363 can be bounded as, The first term in the RHS of Eq (46) is independent of dt and thus is finite as t → ∞. 365 For the second term in RHS of Eq (46), lim t→∞ C 1 (dt) 0 C2(dt) ξ y G(ξ y )dξ y < ∞, but 366 2kBT νdt → ∞ as dt → 0. Since the LHS of Eq (46) is finite, we must have T 1 (R 2 ) = 0.

367
Proof of Theorem 2 We transform Eq (27) into polar coordinates, where the elliptic differential operator and the inner product is defined as, The expression of F r,θ can be found in [26] and we ignore the expression of S for brevity. 373 Then, it's useful to derive F * r,θ (see [26]), the adjoint operator of F r,θ that satisfies 374 v 1 , F r,θ v 2 = F * r,θ v 1 , v 2 , ∀v 1 ∈ D(F r,θ ), v 2 ∈ D(F * r,θ ).