A computational model for the cooperative binding and self-organization of Myosin-VI in membrane reshaping

engineering physics Bachelor’s thesis A computational model for the cooperative binding and self-organization of Myosin-VI in membrane reshaping. Pau Blanco Dorca Supervisor: Prof. Dr. Erwin Frey (LMU Munich) Co-director: Dr. Sergio Alonso Muñoz (UPC) Advisor: Laeschkir Würthner (LMU Munich) June 2020 Acknowledgements I would like to thank Laeschkir Würthner for this guidance throughout this work, also Prof. Dr. Erwin Frey for accepting me to the Statistical and Biological Physics Group in Munich. And I would also like to thank all the members of the group for the welcoming attitude during my brief stay. i Abstract The dynamic membrane reshaping processes are vital for the cells to preform its activity normally. They are involved in a myriad of cellular pathways and can be preformed by different agents. One of those is the myosin-VI, one of the motor proteins that produce the force that will move the membrane. To do that the protein binds to the membrane, this binding has been seen to be targeted towards specific regions of the membrane. In this project a model for the cooperative binding of the myosin VI and a membrane pore will be developed. The growth of the pore will produce a partial derivative equation with moving boundaries that will be solved numerically utilizing the Boundary element method. Then the model will be validated by comparing it to experimental data. ii Contents 1 Introduction 2 1.1 Objectives and work done . . . . . . . . . . . . . . . . . . . . . . . 3 2 Cell biology background 5 2.1 Cellmembrane ............................. 5 2.2 Membrane curvature reshaping . . . . . . . . . . . . . . . . . . . . 6 2.3 Motorproteins ............................. 6 2.3.1 Myosin-VI............................ 8 3 Experimental setup 11 4 The model 14 4.1 Linear stability analysis . . . . . . . . . . . . . . . . . . . . . . . . 17 5 Mathematical derivation of the numerical solution 23 5.1 The Laplace Equation . . . . . . . . . . . . . . . . . . . . . . . . . 23 5.2 Reciprocal relation . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 5.3 Kernelfunction ............................. 25 5.4 Boundary integral solution . . . . . . . . . . . . . . . . . . . . . . . 26 5.5 Boundary solution with constant elements . . . . . . . . . . . . . . 27 5.6 Evolution of the curve . . . . . . . . . . . . . . . . . . . . . . . . . 30 6 Computer implementation 32 7 Results and discussion 35 7.1 Concentration dependency . . . . . . . . . . . . . . . . . . . . . . . 35 7.2 Poregrowth............................... 36 8 Conclusions and further work 39 9 Appendix: Code 41 9.1 Matlab.................................. 41 9.2 C++................................... 47 1 1 Introduction Molecular and cell biology are the perfect targets for the application of computational models. In these fields, the infinity of processes and interactions between the different biomolecules of the cell are of a complexity that makes them impossible to model with a exhaustive representation of all that is happening. Every broken hydrogen bond or the way every protein is folded to preform a certain action would require an astronomical amount of time to compute, so seldom the solution is to use mathematical models that unite millions of different molecular processes to a single variable. This allows for the modelization of complex processes without losing the important features about them. Another degree of complexity is added when we take into account the irregular geometry where all this cellular processes take place. Even the more simple physical phenomena can present emergent behaviours when they occur in a geometrically complex space. This makes mathematical models, that could be easily solved by hand in a regular or highly symmetric environment impossible to solve analytically. In order to obtain results we must turn to the numerical methods and the computer simulations. These methods will be able to solve any model that can be discretized for a computer to do the calculations for us. For this work we wanted a biological problem that would benefit from a computational approach in order to solve it. This is the reason why the growth of a membrane pore triggered by cooperative the binding of the Myosine-VI was chosen. The pores that emerge from this interaction between the membrane lipids and the motor protein do not have a regular and symmetric shape, they instead present irregular shapes that are impossible to calculate analytically. With the huge advances in the computational numeric simulations along the years a great amount of software has been presented to solve all of the typical mathematical problems that can be encountered. This ranges from the simulation of all kinds of fluids to the computation of the structural integrity of solids of a multitude of shapes, among many others. Although these programs are a really useful toolkit, one of the objectives of this thesis was to create the code for the partial differential equation solver by ourselves in order to gain a greater insight in the way it works, 2 both analytically and numerically. Normally for a simple two dimensional partial differential equation can be solved using any commercial software without much more thought. And this is why the model of a growing pore in a cellular membrane was chosen. This problem has the peculiarity that the shape of it’s boundary evolves with time. This means that the simulation had to be build from the start to be able to freely move the boundary and a way to move this boundary needed to be included. 1.1 Objectives and work done The first objective was to develop the model, based on the results from the experimental work of [6]. This would let us see which are be the mathematical tools that will be needed. To do that, different models for diffusion limited aggregation and the solidification where studied[8] and [7], due to its geometrical similarity with the myosin binding process. Even though they belong to very different areas the growth of the interfaces observed in the laboratory resembled the processes depicted in the studied papers. With this base a model was devised. The parameters of the obtained model needed to be fitted to the experimental results and in order to do that we started deriving the linear stability analysis. The results from this analysis helped us choose the values that would lead to the instabilities that appeared in the experiments. Then we discretized the equations from the model using the boundary element method. This would allow us to solve the partial differential equations present in the model. Other numerical methods to solve the other parts of the problem, such as integrals or derivations of different variables that appear in the model. Such as finite differences or Euler integration. Even though these methods are amongst the simpler ones for the integration and derivation of sampled functions they are fast and since they are to be preformed several times every time step the velocity of the method was an important factor. Once the mathematical derivations where done the model needed to be implemented in a computer in order to solve it. Since this simulation would require an important amount of operations the core language had to be a compiled one 3 because of its speed compared to interpreted ones such as Python or Matlab. The chosen language was C++ because of its versatility and the libraries that could be used to solve some parts of the model, such as the linear systems created by the boundary element method and the option to use parallelism to make the project faster. Once the code was working we extracted the results. Because of the numeric demand of this simulation we obtained the variations of the pore growth with he concentration with a simple circular pore. And then we preformed various simulations to determine the different parameters that generated the desired pores. 4 2 Cell biology background 2.1 Cell membrane Lipids or fatty acids are composed of a hydrophilic head and an hydrophobic tail. This asymmetry leads to the auto-organization of groups of lipids when they are in an aqueous medium. The cell membrane is a lipidical bilayer, in other words, a double sheet of fatty acids arranged in order for the hydrophobic tail end to be in between the layers and the hydrophilic head pointing outside. The different lipids are united between them by non-covalent bonds, much weaker than the bonds that form each biomolecule. These weaker bonds allow the components of the membrane to laterally diffuse along the membrane leaflet, which makes membranes a two-dimensional fluid. This diffusivity depends on the composition of the lipids of the membrane and is different for different parts of the membrane or for the membranes of the different organelles. Membranes also contain a myriad of proteins that preform different functions. These proteins control the rate at which certain molecules go through the bilayer and provide binding sites for other proteins to attach to the membrane, among other functions.[cell] Figure 1: Drawing of the membrane, formed by two layers of lipids and the membrane proteins, in red This structure is not exclusive to the cell boundary, it can also be found in different organelles of the cell such as the Golgi apparatus or the mitochondria among others. These organelles have different shapes to fulfill each purpose, these diverse 5 morphologies need different membrane curvatures, and these curvatures have to be dynamically changed in order to carry on the different cellular processes in each organelle. 2.2 Membrane curvature reshaping For pathways such as cell motility and vesicle trafficking the membrane of the cell and some of its organelles must be reshaped dynamically as the process advances. This is a complex process that is carried on by the joint work of different proteins, both embedded in the membrane and exterior ones, lipids and other agents such as actin filaments and other elements of the cytoskeleton. These biomolecules together create a machinery that can sense and reshape the curvature of the membrane to preform the required tasks. This dynamic reshaping can be generated by different agents that work together to obtain different degrees and orientations of membrane curvature. For example by changing the composition of certain membrane lipids bigger heads can be obtained or the expansion of the acyl tails can induce a curvature,[curvature]. Transmembrane proteins such as potassium channels can also act as a wedge that favours certain curvatures because of its asymmetry. This makes the leaflet with the thinner part of the protein a bit smaller and induces a curvature. The cytoskeleton and motor proteins are also primordial for this process. The versatility of the actin cytoskeleton is key to the creation of high curvature regions. This process, like some of the previous ones, requires mechanical force to take place. This force is provided by the motor proteins that move along the citoskeletal fibers and can bind to the membrane. 2.3 Motor proteins Motor proteins differ from other proteins because they can convert the energy obtained from the ATP hydrolysis into movement along a polymer, such as actin or tubulin. There are three different superfamilies of motor proteins, the kinesins and dyneins that move along tubulin filaments, and myosin molecules that additionally move along actin[3]. 6 us a template to compare the results of the model and fit the parameters that cannot be directly extracted from the experiments. We can look for the distance between hotspots or the evolution of the perimeter and the area of the interface to see if the model accurately represents this process. Once we are certain of that we can change the parameters of the model, such as the bulk diffusion or the concentration of myosin molecules to see how these changes affect the growth of these membrane pores. 13 4 The model After studying the results from the experiment we start to design the model, the important features that this needs to include are the cooperative effect that the curvature of the membrane has on the binding of the myosin and the hotspot distance, that is maintained over the growth once a certain radius is reached. In order to maintain this distance in a growing perimeter we will need peaksplitting instabilities to arise as the radius increases. This model is based on diffusion limited solidification models like [8] or[7]. In order to represent the attachment of the myosin-VI to the interface and the evolution of the interface itself observed in the experimental setup we start by defining the interface of the pore as a closed one dimensional curve C,in R2. Then we will define the concentration of myosin in the curve. As seen in the experiments and posed in the paper[6], this concentration will depend on the curvature of the interface. A greater outwards curvature will imply that more myosin will be able to bind to the interface, as shown in Fig. 4. This dependency will be given by the line tension coefficient τthat will establish how strongly the curvature, κ(q) in the equations into the concentration at the interface. c(q) = c0(1 + τκ(q)) ∀q∈C(1) The value c0will be a fraction of the bulk concentration c∞and its value will be determined with the linear stability analysis. In the bulk the myosin will be moving freely in a solution. To simulate this we will assume that the concentration of myosin will diffuse with a coefficient D. Then, since we are only interested in the part that is in contact with the membrane we will project this three dimensional bulk into a two dimensional surface over the lipid bilayer. This diffusion coefficient will be given by the solution used to emulate the cytosol in the experiment. This diffusion has to be defined in a domain, R∈R2 that contents the exterior of C, and will be reigned by the expression ∂tc(q) = D∇2c(q)∀q∈R(2) 14 Then we need to include the effect of force that the myosin applies to the memFigure 6: Radial concentration profile(b.) created by a circular pore (a.) with the depletion zone near the interface at r0and the evolution of the concentration profile for a time differential brane together with the redistribution of the proteins along the pore interface. In order to obtain this equation we will look at a simple pore, a circular one. We will have the concentration profile shown in Fig. 6, with a depletion zone near the membrane pore due to the binding of the cytosol myosin to the pore interface, leading to the minimum concentration of c0right before the interface. The parameter citakes into account the proteins integrated itno the interface. The number of particles must be conserved, so the number of particles attaching to the interface will determine automatically its growth of the area covered by the interface. Given by the next equation D∂rc|r0= (ci−c0)∂tA(3) Then we use that for a circumference the area differential is defined as dA = 2πrdr. Then,we define the normal velocity of the surface as the derivative of respect of the radius. Then, for a circumference of radius r0, Fig. 6 it will be vn=∂tr0. Finally we assume that the concentration in the membrane interface is much larger than the one in the bulk near it, cic0. This is reasonable because of the high cooperativity that the membrane presents. We then obtain the equation for the normal velocity of the interface. vn=D ci ∂nc|r0(4) 15 With this final condition we would have a sufficient model to describe the results of the experiment. Which put together reads c(q) = c0(1 + τκ(q)) ∀q∈C(5) ∂tc(q) = D∇2c(q)∀q∈R(6) c(q) = c∞|q|→∞ (7) vn(q) = D ci ∂nc(q)∀q∈C(8) This model can be further simplified by observing the parameters of the problem. We compute the Péclet number for this equations, this dimensionless number relates the advection and diffusion terms. To do this we compare the timescale of diffusion tdiff ≈r2 0 D(9) The advective term of this problem will be the velocity at which the interface moves, and its timescale will be tgrowth ≈r0 vn (10) Then the dimensionless Péclet number will be the rate between the two timescales Per0=vnr0 D(11) If we input typical values to the parameters obtained from the experiments[6] and bibliography. Such as r0= 1µm,vn∼10−3µm sand D= 1µm2 swe obtain a value of Pe ∼10−3. This means that the diffusion several orders of magnitude faster than the pore growth. Then we will assume that at every time step of the simulation the concentration on the cytosol has already reached the equilibrium position before we evolve the interface. From this result we can make the approximation from [2] to a Laplace equation for the three dimensional bulk projected into the two dimensional plane where the membrane is. 0 = ∇2c(q)∀q∈R(12) 16 With this simplification we will now preform a linear stability analysis to the system in order to see which combinations of parameters lead to a band of unstable modes that will start the growth of the instabilities in the pore interface. 4.1 Linear stability analysis In order to study the stability of our system we will slightly perturb a circular interface with different modes to see its linear evolution[5]. First we search for the different modes that will fit this perturbation. We start with the 2D Laplace equation in polar coordinates 1 r∂r[r∂rf(r, θ)] + 1 r2∂2 θf= 0 (13) Then we assume that the solution can be obtained by separation of variables. So the function we are looking for will be of the form f(r, θ) = R(r)Θ(θ). If we put this function into the Laplace equation in polar coordinates we obtain Θ∂2 rR+ Θ∂rR r+∂2 θΘR r2= 0 (14) By reordering the equation we obtain r2∂2 rR R+r∂rR R=−∂2 θΘ Θ(15) Since this equality holds for any point we can assume that both parts are constant, so we define the constant n, that will determine the mode we are studying r2∂2 rR R+r∂rR R=n2∂2 θΘ Θ=−n2(16) Both differential equations can be easily solved and we obtain R(r) = r−nΘ(θ) = cos nθ (17) Which are not the unique solutions but are enough for the problem at hand, but they are enough to describe the modes that we will use. If we join both the radial 17 and the angular parts we obtain a series of solutions that will be named Pn Pn(r, θ) = r−ncos nθ (18) Then we can define the perturbed interface as a circle with a modulation defined by the mode r(θ) = r0+δcos nθ (19) With r0the radius of the initial circular interface and δthe amplitude of the small perturbation, it is assumed that r0δ. Once we have defined the interface we can calculate the boundary condition as stated in the model Eq. (1). For this we need to calculate the curvature of the interface, which we can obtain using the following formula, where the dot notation refers to the derivative with respect to θ κ(θ) = |r2+ 2 ˙r2−r¨r| (r2+ ˙r2)3/2(20) Since we are working with a slightly perturbed sphere we will consider the terms of the derivative squared to be negligible against the squared value of the radius and that the second derivative of the radius will be always smaller than the radius itself. With this approximations we obtain κ(θ)≈r2−r¨r r3=1 r−¨r r2(21) Then we can approximate the first term around r0to the first order of δ 1 r=1 r0 −1 r2 0 (r−r0) = 1 r01−δcos nθ r0+O(δ2)(22) Then we can compute ¨rby using Eq. (19) ¨r=δcos nθ =−n2δcos nθ (23) 18 With this result we can also approximate the second term −n2δcos nθ r2≈ −n2δcos nθ 1 r2 0 −2 r3 0 (r−r0)=−n2δcos nθ r2 01−2δcos nθ r0≈ −n2δcos nθ r2 0 (24) Once we have done all the approximations to the first order of δwe obtain the curvature κ(θ)≈1 r01−δcos nθ r0+n2δcos nθ r2 0 =1 r01 + δcos nθ r0n2−1(25) Then by using this result and Eq. (1) the concentration in the interface Cwill be cC(θ)≈c01 + τ r01 + δcos nθ r0 (n2−1) (26) Then we can compute the solution on the rest of the domain R. In order to do this we have to provide a solution of the Laplace equation that will be equal to Eq. ( 26) at the interface and also fulfill the condition c(r1) = c∞with r1r0. This solution will have the general form of c(r, θ) = Aln r+BδPn(r, θ) + C(27) By expanding it to the first order of δaround the interface we obtain cC(θ) = Aln r0+δcos nθ r0+Bδ cos nθ rn 0 +C(28) If we equal this last equation to Eq. (26) and use the boundary condition for the bulk we obtain the three coefficients A=1 lG+c0τ r0B=rn−1 0mc0τ r0n2−1−1 l−G lC=c∞−Aln r1 (29) With G=c0−c∞l= ln r0 r1 (30) Once we have the concentration field we have to compute its normal derivative in 19 the interface to see the evolution of its perimeter. Since the interface has an almost circular shape we will approximate the normal derivative as a radial derivative ∂rc. Then by taking the the radial derivative of Eq. (27) we obtain ∂rc=A r−Bnδ cos nθ rn+1 (31) And by approximating it around r0we obtain ∂rc≈1 lr0G+c0τ r0−δcos nθ r2 0nc0τ r0n2−1+G+c0τ r01−n l(32) Then by using the definition of the normal velocity of the interface given in the model we obtain vn=D ci1 lr0G+c0τ r0−δcos nθ r2 0nc0τ r0n2−1+G+c0τ r01−n l (33) This result gives us the velocity of the whole interface, but we can extract the velocity of a mode by taking the time derivative of Eq. (19) vn=dr(θ) dt =dr0 dt +dδ dt cos nθ (34) This result lets us separate the growth of the unperturbed circle and the growth of the perturbation. The radius of the circle will grow following the next expression ˙r0=D lcir0G+c0τ r0(35) This result yields that the initial interface needs to have a minimum radius in order to start its growth, this critical radius rcis given by rc=c0τ c1−c0(36) And the evolution of the perturbation amplitude will be given by the next expression ˙ δ δ=−D cir2 0nc0τ r0n2−1+G+c0τ r01−n l(37) 20 Once we know the linear behaviour we can obtain the dispersion relation for different groups of parameters, this will help us chose the correct parameters for the simulations. We start by setting some of the known parameters from the experiments. We will set the concentration from a range around 50nM to 150nM. Then we will set the initial radius to 300nm to emulate the initial dimensions of the pores. With the bulk concentration and the radius fixed, at least within a range, it is observed that only the ratio between c0and c∞, and the value of the line tension τwill change the stability of the system. Furthermore approximate threshold values for these parameters can obtained from the hotspot distance as a function of the diameter, in Fig. (5). Then by assuming a circle-like shape we can establish that the highest unstable mode must fulfill nmax >πd hdzi(38) With hdzithe hotspot distance. By looking at the figure we can obtain the minimum value for the highest mode that will aid with the posterior choice of parameters r0(µm)nmax 0.5 8 1 11 1.5 14 2 16 2.5 18 3 19 3.5 22 4 25 These values are approximated and are only valid for the lineal regime, but they serve as a guideline to choose the range of parameters that we will focus on. We can see from the results that the fist mode will always maintain the same amplitude, this is coherent with the fact that it only corresponds to a translation of the whole interface, and the center of interface will remain in the same place for any sinusoidal mode. It can also be observed that the increase of the radius will 21 Figure 7: Dispersion relation for different values of τmaintaining the rest of parameters constant. It can clearly be seen that if we increase this parameter the number of unstable frequencies diminishes increase the unstable modes. This is coherent with the peak splitting instabilities that we want to obtain, since this higher modes will appear as the radius grows. 22 By using this approximations we can rewrite [57] as c(p) = N X k=1 [ckIk 2(p)−vkIk 1(p)] ∀p∈R(62) 1 2c(p) = N X k=1 [ckIk 2(p)−vkIk 1(p)] ∀p∈R(63) Ik 1(p) = ZCk Φ(q, p)dsq(64) Ik 2(p) = ZCk ∂nkΦ(q, p)dsq(65) Then we can determine the values of the normal derivative of the function at the boundary by evaluating [62] at the midpoints of the segments. 1 2cm= N X k=1 [ckIk 2(dm)−vkIk 1(dm)] (66) In order to obtain the values we will take m= 1, ..., N and obtain the following linear system M1·v=M2·c(67) M1m k=Ik 1(dm)(68) M2m k=Ik 2(dm)−1 2δm k(69) With δm kKronecker delta. Then if we solve for vwe obtain the normal derivative at the boundary that will allow us to obtain the normal velocity of the interface v=M1 −1M2c(70) Once we know both ckand vkwe can use Eq. (62) with λ= 1 to compute the values ∀p∈R c(p)≈ N X k=1 [ckIk 2(p)−vkIk 1(p)] (71) 29 5.6 Evolution of the curve Variables are labeled with a subindex iif they are evaluated at the beginning of the segment, at pi, and a subindex jif they are evaluated in the middle point, di. To convert from a function evaluated at the beginning of the segments fito a function evaluated at the midpoints fjwe use fj=fi+fi+1 2fi=fj−1+fj 2(72) With this nomenclature, the variables that we obtain from the boundary element method are ct jand vt j. This values provide enough information evolve the curve for one time step. The evolution of the curve is given only by the normal velocity of the curve, as exposed in [6]. To evolve the curve we will use the the integral method from the aforementioned paper. We can compute the evolution of the perimeter by using the following integral ∂tL=Z1 0 vn∂ρθdρ (73) With ρan arc lent parametrization escalated to range from 0 to 1 along the curve. Now we need to define a discrete integration method. Since the evaluation of the functions is in the middle points the numerical integration will be trapezoidal Then, any integral will become Zb a f(t)dt ≈Intt[fi]b a(74) With Intt[fi]b a= N−1 X i=0 (fi+fi+1)(ti+1 −ti) 2= N−1 X i=0 fj(ti+1 −ti) 2(75) For a vector with Nsamples. If we use the results obtained from solving the PDE the discretized derivative of the perimeter becomes dtLt=Intρ[vt i∂ρθt i]1 0(76) 30 Then Lt, the discretized version of L(t), is obtained by using a simple Euler method Lt=Lt−∆t+ ∆t·dtLt−∆t(77) Once we have the updated length we can compute the evolution of the angle of the segments of the interface with the positive horizontal axis θ(ρ). First we obtain the partial derivative ∂tθ(t, ρ)given by ∂tθ(t, ρ) = −1 L(t)∂ρVn(t, ρ) −1 L(t)∂ρθρZ1 0 Vn(t, ρ)∂ρθdρ −Zρ 0 Vn(t, ρ)∂ρθdρ(78) In the discretized notation will become dtθt i=−1 Ltdρvt ni −1 LtdρθiρiIntρ[vt nidρθt i]1 0−Intρ[vt nidρθt i]ρi 0(79) Then we evolve the discretized θt iusing again the Euler integration. θt i=θt−∆t i+ ∆tdtθt−∆t i(80) Once we have the updated theta we can retrieve the values of the positions of the interface by using the following integrals x(t, ρ) = x(0,0) + Zt 0 vn(t0,0) sin θ(t0,0)dt0+L(t)Zρ 0 cos θ(t, ρ0)dρ0(81) y(t, ρ) = y(0,0) −Zt 0 vn(t0,0) cos θ(t0,0)dt0+L(t)Zρ 0 sin θ(t, ρ0)dρ0(82) Which will be discretized in the following way xt i=x0 0+Intt[vt0 n0sin θt0 0]t 0+LtIntρ[cos θt i]ρ 0(83) yt i=y0 0−Intt[vt0 n0cos θt0 0]t 0+LtIntρ[sin θt i]ρ 0(84) This leads to the new position of the interface which can be used to solve again the Laplace equation and repeat the process for the next time step. 31 6 Computer implementation With the discretized model fully developed we start to implement it to a computer. To do that we need to choose a computer language to convert it. After some consideration C++ was chosen for the computation of the evolution of the interface for its speed and versatility. From the mathematical derivation it is clear that a library for the creation and arithmetics will be needed. In order to supply that, the Eigen library was used. With this library the creation and algebraic manipulation of matrices were much more straightforward. Since we are working with an infinite two-dimensional bulk Figure 10: Points of the discretitzed boundary where the concentration is known. This includes both the interface boundary and points that set the concentration of the solution far from the interface. Representation of the space that these points fill in the Boundary element matrices, this matrices compute the influence of each point over all the others. Only the yellow part is constant. The blue and mixed color parts will change as the interface evolves since the distance between the points will change, meaning that the influence they exert to each other will also be modified. 32 we need to bound it at a particular radius in order to set the points where we set the concentration of the bulk. This points will be set in along a circumference of radius a hundred times bigger than r0, as seen in the linear stability analysis. Then the collocation points will be as shown in Fig [10]. To solve the system of linear equations created by the boundary element method, we chose the Householder QR decomposition algorithm, provided by the Eigen library, which offered the best trade-off between velocity and accuracy. And a fast matrix solver is an important aspect for this particular problem. This is because for our moving boundary problem almost all of the matrix needed to be recalculated at every time step, as seen in Fig. [10]. This means that if the matrix solving algorithm is not efficient the time of the simulation would be impractically high and the number of points used to discretize the boundary would need to be reduced. This carries with it problems since the distance between the points of the boundary has to be small in order to avoid numeric instabilities. This problem also gets aggravated by the fact that the number of points to which the interface is discretized also grows with the pore perimeter. This is done to maintain the distance between consecutive points bounded within the numerically stable range. To do this every time the distance between two consecutive points is greater than a fixed threshold we interpolate the the known positions of the interface to double the number of points. Taking this limitations into account we wrote the code that would solve and evolve the model. And to test it we started by trying the growth of single mode interfaces. In order to keep all the points equispaced in the interface and minimize the numerical instabilities we set the initial condition analogously to [[6]] and set the initial condition for the angle of each point to θn(ρ)=2πρ +δsin(2πnρ)(85) With nthe mode and ρthe parametrization of the curve, ranging from zero to one. Then scale the points to the initial radius and run and start advancing time. This will help choose the remaining parameters in order to obtain a growth similar to the one seen in the experiments. 33 Once the model is tested we set different initial conditions, we started with a sum of modes to see the competition between them. θ(ρ) = 2πρ + N X n=1 δ(n) sin(2πnρ +u(n)) (86) Where we create an initial interface that contains a sum of the first the first N modes, with a random phase induced by the uniform distribution u(n), that goes from 0to π. These modes will vanish or increase and interact between them, giving . In the next section we will review the results of these simulations. Finally we also used a circle perturbed by a random Gaussian noise that would better represent the irregularities in the membrane that would act as seed for the growth of the pore. in the form of θ(ρ) = 2πρ +gδ(ρ)(87) With gδa random number extracted from a Gaussian distribution centered at zero and with a standard deviation of δ. By repetition of the simulations and analysis and comparison of the results with the experiments a set of parameters was chosen. Some of them were given by the experimental results and the other were obtained from fitting the results of the simulations. 34 7 Results and discussion Once the parameters of the model are chosen we can start to extract results from the simulations. These will be the variation in the growth with the concentration and the evolution of the pore to, taking special attention to the splitting of the hotspots with the growth of the radius due to the growth of the unstable region seen in the linear stability analysis. 7.1 Concentration dependency Figure 11: Evolution of the perimeter for different concentrations, it can be seen that the velocity of the interfaces grows with the concentration without changing its profile One of the effects that couldn’t be obtained from the phenomenological model in [6] was the dependency of the growth of the pores with the concentration. The experimental data clearly depicted an increase of the velocity of growth with the increase of myosin VI in the cytosol. One of the main objectives of this work was to obtain a dependency with the concentration that resembled the one seen in the experiments. To do this we created circular pores, which needed less computational time to grow since they do not need many points for the solutions to be numerically stable. 35 This yields the results shown in Fig. (11), that can be compared to the experimental ones in Fig. (5). We observe that the changes in concentration create a change in the timescale of the whole problem, augmenting the velocity of the process for higher concentrations. Through this comparison the parameters such as the relation between c0and c∞, and cicould be further polished to better fit the experimental results available. 7.2 Pore growth Since this is a non-local model the whole interface contributes to the evolution of every point of the interface. This hinders the fitting of the parameters since for different initial seeds the evolution of the interface can vary significantly. In order to obtain these parameters, the ones that are not fixed, like the concentration or the initial radius, several initial conditions have been tested. The parameters for each of those where chosen according to the linear stability analysis preformed. Each seed was initialized with a slightly different set of parameters, shown below,that allowed us to observe different behaviours. Seed number 1 2 3 4 5 6 c00.5 0.5 1 1 0.5 0.5 c∞50 50 50 50 50 50 D c017 17 17 17 17 17 r0200 200 200 300 200 200 τ10 4 10 5 4 4.5 Then these results where compared with the with the experiments in order to see if they could accurately model them. It can be seen, Fig (13), that the radius versus perimeter power law is fulfilled in all the cases. The distance between hotspots fluctuates more aggressively than the experimental data. This may be because the exactly correct seed to start the simulation hasn’t been found. On the other hand is clear that seeds with large curvature changes, like the second and the sixth create a higher first spike. Nevertheless the overall trend seems to be stabilizing at a fixed value even with the initial instabilities. In this plot the fifth seed wasn’t 36 Figure 12: Results for the evolution of the pore. In the figure it can be seen the importance of the initial conditions in the growth of the pore. The seeds with higher mode perturbations such as 4. will led to more dramatic inward curvature. And the seeds with lesser modes, such as 5. will no provide instabilities to accurately model growth of the pore. 2. and 6. have the same initial condition but the increased τin 6 will make it harder for the instabilities to appear. Finally we have 1. the sum of different modes and 3. with a seed that is perturbed by a random Gaussian noise. The irregularity of the noise breaks the symmetry seen in the other flowers. included. Finally for the concentration we can see that in general the perimeter of the pores grows at the expected rate. The only one that seems to deviate more from this trend is the fifth seed, whose oval shape kept the the perimeter from growing by avoiding the inward curvature and peak splitting present in the other pores. 37 Figure 13: Results of the simulation for the six different initial conditions compared against the experimental data. The plots of the experimental data have been cleaned to avoid clumping. To see the whole information refer to Fig. (5). 38 F2(mm,kk) = 0; else ff1 = B/A; ff2 = E/A; ff3 = sqrt(abs(ff3)); f21 = l(kk)*(n(kk,1)*(p(kk,1)-cm(mm,1)) ... + n(kk,2)*(p(kk,2)-cm(mm,2)))/(pi*ff3); f22 = atan((2*A+B)/ff3)-atan(B/ff3); F2(mm,kk) = f21*f22; F1(mm,kk) = l(kk)/(4*pi)*(2*(log(l(kk))-1)- ...ff1/2*log(abs(ff2))+(1+ff1/2) ...*log(abs(1+ff1+ff2))+ff3/A*f22); end end end end % BEM solve for mm = 1:N for kk = 1:N M1(mm,kk) = F1(mm,kk); M2(mm,kk) = F2(mm,kk) - 0.5*eq(mm,kk); end end b = M2*bc; c = M1\b; for ii = 1:N phi(ii) = bc(ii); v(ii) = c(ii); end % Vel dL = zeros(N1+1,1); 45 vnj = -aa^2*d*v(1:N1); drhoj = l(1:N1)/L; dkj = (circshift((dthj/L),-1) + circshift((dthj/L),1) ...-2*(dthj/L))./(l(1:N1).^2); dkj = 0.25*(2*dkj + circshift(dkj,1) + circshift(dkj,-1)); vnj = vnj + k1*dkj; vni = 0.5.*(vnj + circshift(vnj,1)); dvnj = (circshift(vni,-1)-vni)./drhoj(1:N1); rhoi = [0;cumsum(l(1:N1))/L]; rhoj = 0.5.*(rhoi + circshift(rhoi,-1)); rhoj = rhoj(1:N1); dL(1) = 0; for ii = 2:N1+1 dL(ii) = dL(ii-1) + drhoj(ii-1)*dthj(ii-1)*vnj(ii-1); end dLj = 0.5.*(dL + circshift(dL,-1)); dLj = dLj(1:N1); ff = rhoj.*dL(end) - dLj; dtthj = 1/L*dvnj -1/L*dthj.*ff; thj = thj + dt*circshift(dtthj,0); thi = 0.5*(thj + circshift(thj,1)); thi(1) = thi(1) + pi; ll(tt) = L; pn(1,:) = [(p(1,1) - dt*vni(1)*sin(thi(1))) ...(p(1,2) + dt*vni(1)*cos(thi(1)))]; L = L + dt*dL(end); for ii = 2:N1 pn(ii,:) = pn(ii-1,:) + L*[cos(thj(ii-1))*drhoj(ii-1) ...sin(thj(ii-1))*drhoj(ii-1)]; end 46 p(1:N1,:) = pn; if mod(tt,1000) == 0 stg = stg +1; stgthi{stg} = 0.5*(dthj + circshift(dthj,-1)); stgp2{stg} = pn; end tt = tt + 1; ll(tt) = L; init = 0; end end } 9.2 C++ #include <iostream> #include <fstream> #include <Eigen/Dense> #include <math.h> #include <random> using namespace std; using Eigen::MatrixXd; using Eigen::ArrayXXd; using Eigen::ArrayXd; using Eigen::VectorXd; ////////////////// // Declarations // ////////////////// // Parameter definition 47 const int N1 = 500; // Number of elements in the interface const int N2 = 50; // Number of elements in the bulk boundary const int N = N1 + N2; // Total number of points const long double r1 = 216; // Initial radius of the interface const long double r2 = r1*100; // Radius of the bulk points const long double tau = -1; // Line tension const long double c0 = 0.1; // Bulk concentration int ttf = 10000; // Steps of the simulation long double dt = 0.0001; // Time step const long double tf = dt*ttf; // Time of simulation const long double stddev = 0.0; long double dp = 0.05; //Amplitude of the perturbation long double Ln,L,dL,ff; int ttstg = 0; // Number of times that the perimeter is saved int tt = 0; // Current iteration number int tt0 = 0; // bool init = 1;// If saved initial conditions init = 0, else init = 1 long double xx0, yy0; const long double a = 10,d = -1; double rnd,md; int mdmax = 20; double xin,yin; default_random_engine gen; normal_distribution<long double> noise(r1,stddev); // Variable declaration ArrayXXd p(N,2), n(N,2),cm(N,2); ArrayXd lj(N),vnj(N1), dvnj(N1), vni(N1), rhoj(N1), ...rhoi(N1), thj(N1), dthi(N1), ...dthj(N1), unif(ttf+N); 48 ArrayXd dtthj(N1), dtthi(N1), thi(N1), thin(N1), mmj(N1); ArrayXd ll(ttf+1),vnit(ttf),thit(ttf), mdr(mdmax); MatrixXd F1(N,N), F2(N,N); VectorXd cj(N), dcj(N), mm(N); /////////////// // Functions // /////////////// // Basic integration function double integ(int N, ArrayXd lj,ArrayXd fi) { double in = 0; for (int ii = 0; ii < N-1; ii++) { in += 0.5*lj(ii)*(fi(ii) + fi(ii+1)); } return in; } double integj(int N, ArrayXd lj,ArrayXd fj) { double in = 0; for (int ii = 0; ii < N; ii++) { in += lj(ii)*fj(ii); } return in; } // Basic cosine integration function double cinteg(int N, ArrayXd lj,ArrayXd fi) { 49 double in = 0; for (int ii = 0; ii < N-1; ii++) { in += lj(ii)*cos((fi(ii) + fi(ii+1))/2); } return in; } // Basic cosine integration function double sinteg(int N, ArrayXd lj,ArrayXd fi) { double in = 0; for (int ii = 0; ii < N-1; ii++) { in += lj(ii)*sin((fi(ii) + fi(ii+1))/2); } return in; } // From point to midpoint ArrayXd itoj(int N, ArrayXd fi) { ArrayXd fj(N); for (int ii = 0; ii < N-1; ii++) { fj(ii) = 0.5*(fi(ii) + fi(ii+1)); } fj(N-1) = 0.5*(fi(0) + fi(N-1)); return fj; } // From midpoint to point ArrayXd jtoi(int N, ArrayXd fj) { ArrayXd fi(N); 50 for (int ii = 1; ii < N; ii++) { fi(ii) = 0.5*(fj(ii) + fj(ii-1)); } fi(0) = 0.5*(fj(0) + fj(N-1)); return fi; } // Unwraping function for the angular variables ArrayXd unwrap(int N, ArrayXd f) { bool wrap = 1; while (f(0) > 6) { f(0) -= 2*M_PI; } while (f(0) < 3) { f(0) += 2*M_PI; } while (wrap == 1) { wrap = 0; for (int ii = 1; ii < N; ii++) { if (f(ii)-f(ii-1) > 4) { f(ii) -= 2*M_PI; wrap = 1; } else if (f(ii)-f(ii-1) < -4) { f(ii) += 2*M_PI; wrap = 1; } } } return f; } 51 ///////////////// // Main script // ///////////////// int main() { // Openig files for saving data ofstream stgdata; if (init == 1) { stgdata.open("stgdata.txt",ofstream::out); } else { stgdata.open("stgdata.txt",ofstream::out | ofstream::app); } // Initial conditions if (init == 1) { for (int jj = 0; jj < mdmax; jj++){ mdr(jj) = 100*noise(gen); } for (int ii = 1; ii<N1 + 1; ii++) { // Circle + noise rnd = noise(gen); p.row(N1-ii) << rnd*cos(2*ii*M_PI/N1),rnd*sin(2*ii*M_PI/N1); //Circle //p.row(N1-ii) << r1*cos(2*ii*M_PI/N1),r1*sin(2*ii*M_PI/N1); //Circle + mode p.row(N1-ii) << r1*(1+dp*cos(8*ii*M_PI/N1))*cos(2*ii*M_PI/N1),r1* ...(1+dp*cos(8*ii*M_PI/N1))*sin(2*ii*M_PI/N1); 52 //Circle + modes md = 0; for (int jj = 1; jj <mdmax; jj++){ md += dp*cos(2*jj*ii*M_PI/N1); } //p.row(N1-ii) << rnd*(1+md)*cos(2*ii*M_PI/N1),rnd*(1+md)*sin(2*ii*M_PI/N1); } for (int ii = 1; ii<N2 + 1; ii++) { p.row(N1 + ii - 1) << r2*cos(2*ii*M_PI/N2),r2*sin(2*ii*M_PI/N2); } } else { cin >> tt0; for (int ii = 0; ii < N ; ii++) { cin >> xin >> yin; p.row(ii) << xin,yin; } } xx0 = p(0,0); yy0 = p(0,1); for (int ii = 0; ii<ttf+N; ii++) { unif(ii) = 1; } // Time evolution while (tt < ttf) { for (int ii = 0; ii<N1-1; ii++) { lj(ii) = sqrt(pow(p(ii+1,0)-p(ii,0),2) + pow(p(ii+1,1)-p(ii,1),2)); n.row(ii) << (p((ii+1),1)-p(ii,1))/lj(ii),(p(ii,0)-p((ii+1),0))/lj(ii); cm.row(ii) << (p(ii,0) + p(ii+1,0))/2,(p(ii,1) + p(ii+1,1))/2; } lj(N1-1) = sqrt(pow(p(0,0)-p(N1-1,0),2) + pow(p(0,1)-p(N1-1,1),2)); 53 n.row(N1-1) << (p(0,1)-p(N1-1,1))/lj(N1-1),(p(N1-1,0)-p(0,0))/lj(N1-1); cm.row(N1-1) << (p(N1-1,0) + p(0,0))/2,(p(N1-1,1) + p(0,1))/2; for (int ii = N1; ii<N-1; ii++) { lj(ii) = sqrt(pow(p(ii+1,0)-p(ii,0),2) + pow(p(ii+1,1)-p(ii,1),2)); n.row(ii) << (p((ii+1),1)-p(ii,1))/lj(ii),(p(ii,0)-p((ii+1),0))/lj(ii); cm.row(ii) << p.row(ii)/2 + p.row(ii + 1)/2; } lj(N-1) = sqrt(pow(p(N1,0)-p(N-1,0),2) + pow(p(N1,1)-p(N-1,1),2)); n.row(N-1) << (p(N1,1)-p(N-1,1))/lj(N-1),(p(N-1,0)-p(N1,0))/lj(N-1); cm.row(N-1) << p.row(N-1)/2 + p.row(N1)/2; // Curvature of the interface if(tt == 0){ L = integ(N1+1,lj,unif); for (int ii = 0; ii<N1; ii++) { if (p(ii+1,1) - p(ii,1) > 0){ if (p(ii+1,0) - p(ii,0) > 0){ thj(ii) = acos((p(ii+1,0)-p(ii,0))/lj(ii)); } else { thj(ii) = M_PI - acos(-(p(ii+1,0)-p(ii,0))/lj(ii)); } } else { if (p(ii+1,0) - p(ii,0) > 0) { thj(ii) = -acos((p(ii+1,0)-p(ii,0))/lj(ii)); } else { thj(ii) = M_PI + acos(-(p(ii+1,0) - p(ii,0))/lj(ii)); } } } 54