Parametric FEM for Shape Optimization applied to Golgi Stack

TL;DR

This paper develops a parametric finite element method for shape optimization of Golgi stacks, modeling cisternae as elastic vesicles to understand their morphological evolution.

Contribution

It introduces a novel FEM-based shape optimization model for Golgi stacks, incorporating biological constraints and handling multiple vesicles simultaneously.

Findings

The model successfully simulates Golgi stack shapes.

Numerical results align with observed Golgi morphologies.

The approach offers insights into Golgi cisternae dynamics.

Abstract

The thesis is about an application of the shape optimization to the morphological evolution of Golgi stack. Golgi stack consists of multiple layers of cisternae. It is an organelle in the biological cells. Inspired by the Helfrich Model \cite{Helfrich}, which is a model for vesicles typically applied to biological cells, a new model specially designed for Golgi stack is developed and then implemented using FEM in this thesis. In the Golgi model, each cisternae of the Golgi stack is viewed as a closed vesicle without topological changes, and our model is adaptable to both single-vesicle case and multiple-vesicle case. The main idea of the math model is to minimize the elastic energy(bending energy) of the vesicles, with some constraints designed regarding the biological properties of Golgi stack. With these constraints attached to the math model, we could extend this model to an obstacle-type problem. Hence, in the thesis, not only the simulations of Golgi stack are shown, but some interesting examples without biological meanings are also demonstrated. Also, as multiple cisternaes are considered as a whole, this is also a model handling multiple objects. A set of numerical examples is shown to compare with the observed shape of Golgi stack, so we can lay down some possible explanations to the morphological performance of trans-Golgi cisternae.

Figures

Tables

FEM	Finite Element Method
SDG	Staggered Discontinuous Galerkin

$ℝ$	the set of real numbers
$\Gamma$	a hypersurface of dimension $k-1$ embedded in ${ℝ}^k$
$𝝂$	outer unit normal vector field
$\parallel \cdot \parallel$	Euclidean norm
${\mathbb{P}}_2\left(K\right)$	The space of polynomials of degree $\le$ 2 defined over the set $K$
$C^1$	continuously differentiable functions
$D\left(\varphi \right)$	$D{\left(\varphi \right)}_{ij}={\left[{\nabla }_{\Gamma }{\varphi }_j\right]}_i+{\left[{\nabla }_{\Gamma }{\varphi }_i\right]}_j$
${\nabla }_{\Gamma }\varphi$	tangential gradient of $\varphi$
${𝒟}_{\Gamma }\left(\varphi \right)$	tangential Jacobian matrix of $𝐰$
${\text{div}}_{\Gamma }\varphi$	tangential divergence of $𝐰$
${\Delta }_{\Gamma }\varphi$	tangential Laplace of $\varphi$

Equations

$$E(\Gamma )=\frac{K_b}{2}{\int }_{\Gamma }{h}^{2}dS+\frac{K_G}{2}{\int }_{\Gamma }GdS,$$

$${\int }_{\Gamma }GdS=2\pi \chi (\Gamma ),$$

$$W(\Gamma ):={\int }_{\Gamma }{h}^{2}dS,$$

$$A(\Gamma ):={\int }_{\Gamma }1d\Gamma ,$$

$$H(\Gamma ):={\int }_{\Gamma }{1}_{\mathbf{B}}(\mathbf{x})d\Gamma ,$$

$${\nabla }_{\Gamma }\varphi :=\nabla \tilde{\varphi }{|}_{\Gamma }-\frac{\partial \tilde{\varphi }}{\partial n}n,$$

$$[{\mathcal{D}}_{\Gamma }(\mathbf{\varphi }){]}_{ij}=[{\nabla }_{\Gamma }{\mathbf{\varphi }}_i{]}_j.$$

$${\text{div}}_{\Gamma }(\mathbf{\varphi }):=\text{tr}({\mathcal{D}}_{\Gamma }(\mathbf{\varphi })),$$

$${\Delta }_{\Gamma }\varphi ={\text{div}}_{\Gamma }{\nabla }_{\Gamma }\varphi$$

$$J(\Gamma )={\int }_{\Gamma }\varphi d\Gamma .$$

$$dJ(\Gamma ;\mathbf{v})=\underset{t\to 0^{+}}{\lim }\frac{J({\Gamma }_t)-J(\Gamma )}{t},$$

$$\begin{array}[]{ll}dW(\Gamma;\bm{v})=\int_{\Gamma}\nabla_{\Gamma}{\bm{v}}\colon\nabla_{\Gamma}{\bf h}&-\int_{\Gamma}\nabla_{\Gamma}\bm{v}(\nabla_{\Gamma}\text{Id}+\nabla_{\Gamma}\text{Id}^{T})\colon\nabla_{\Gamma}{\bf h}\\ &+\frac{1}{2}\int_{\Gamma}\text{div}_{\Gamma}{\bf h}\ \text{div}_{\Gamma}\bm{v}\end{array}$$

$$\begin{array}[]{ll}dW(\Gamma;\bm{v})=&-\int_{\Gamma}\nabla_{\Gamma}\bm{v}\colon\nabla_{\Gamma}{\bf h}+\int_{\Gamma}D(\bm{v})\nabla_{\Gamma}\text{Id}\colon\nabla_{\Gamma}{\bf h}\\ &-\int_{\Gamma}\text{div}_{\Gamma}{\bf h}\ \text{div}_{\Gamma}\bm{v}-\frac{1}{2}\int_{\Gamma}|h|^{2}\text{div}_{\Gamma}\bm{v}\end{array},$$

$$J(\Gamma )={\int }_{\Gamma }\varphi d\Gamma .$$

$$dJ(\Gamma ;\mathbf{v})={\int }_{\Gamma }\nabla \varphi \cdot \mathbf{v}+\varphi ({\text{div}}_{\Gamma }\mathbf{v})d\Gamma .$$

$$dA(\Gamma ;\mathbf{v})={\int }_{\Gamma }{\text{div}}_{\Gamma }\mathbf{v}d\Gamma .$$

$$dA(\Gamma ;\mathbf{v})={\int }_{\Gamma }\mathbf{v}\cdot \mathbf{h}d\Gamma .$$

$$⟨\mathbf{v},\varphi ⟩=-dJ(\Gamma ;\varphi ),\forall \varphi \in Hil(\Gamma ),$$

$$\begin{array}[]{cc}\displaystyle\text{minimize}&W(\Gamma)=\frac{1}{2}\int_{\Gamma}h^{2}\,d\Gamma,\\ \text{subject to}&A(\Gamma)=A(\Gamma_{0}),\end{array}$$

$$J_1(\lambda ,\Gamma ):=W(\Gamma )+\lambda (A(\Gamma )-A({\Gamma }_0))=\frac{1}{2}{\int }_{\Gamma }{h}^{2}d\Gamma +\lambda \left({\int }_{\Gamma }d\Gamma -{\int }_{{\Gamma }_0}d{\Gamma }_0\right).$$

$$G_T=\bigcup _{t\in [0,T]}\Gamma (t)\times \{t\}$$

$$\mathbf{u}(x,t)=x$$

$${\int }_{\Gamma (t)}\dot{\mathbf{u}}\cdot \varphi =-d{J}_1(\lambda (t),\Gamma (t);\varphi )=-dW(\Gamma (t);\varphi )-\lambda (t)dA(\Gamma (t);\varphi ),$$

$$A(\Gamma (t))=A({\Gamma }_0),$$

$$\begin{array}[]{ll}dW(\Gamma;\phi)=\int_{\Gamma}\nabla_{\Gamma}{\bf\phi}\colon\nabla_{\Gamma}{\bf h}&-\int_{\Gamma}\nabla_{\Gamma}{\bf\phi}(\nabla_{\Gamma}\text{Id}+\nabla_{\Gamma}\text{Id}^{T})\colon\nabla_{\Gamma}{\bf h}\\ &+\frac{1}{2}\int_{\Gamma}\text{div}_{\Gamma}{\bf h}\ \text{div}_{\Gamma}{\bf\phi},\end{array}$$

$$\begin{array}[]{ll}dW(\Gamma;\phi)=&-\int_{\Gamma}\nabla_{\Gamma}{\bf\phi}\colon\nabla_{\Gamma}{\bf h}+\int_{\Gamma}D(\phi)\nabla_{\Gamma}\text{Id}\colon\nabla_{\Gamma}{\bf h}\\ &-\int_{\Gamma}\text{div}_{\Gamma}{\bf h}\ \text{div}_{\Gamma}{\bf\phi}-\frac{1}{2}\int_{\Gamma}|h|^{2}\text{div}_{\Gamma}{\bf\phi},\end{array}$$

$$dA(\Gamma ;\varphi )={\int }_{\Gamma }\mathbf{h}\cdot \varphi .$$

$$H(\Gamma )={\int }_{\Gamma }{1}_{\mathbf{B}}(x)d\Gamma$$

$$\begin{array}[]{cc}\displaystyle\text{minimize}&W(\Gamma)+\alpha H(\Gamma),\\ \text{subject to}&A(\Gamma)=A(\Gamma_{0}).\end{array}$$

$$\begin{array}[]{ll}J_{2}(\lambda,\Gamma)&:=W(\Gamma)+\alpha H(\Gamma)+\lambda(A(\Gamma)-A(\Gamma_{0}))\\ &=\frac{1}{2}\int_{\Gamma}h^{2}\,d\Gamma+\alpha\int_{\Gamma}1_{\bf B}(x)\,d\Gamma+\lambda\left(\int_{\Gamma}\,d\Gamma-\int_{\Gamma_{0}}\,d\Gamma_{0}\right),\end{array}$$

Full text

1. 1 Introduction

2. 1.1 Overall Introduction and Motivation

3. 1.1.1 Willmore Energy

4. 1.1.2 Area Constraint

5. 1.1.3 Barrier Functional

6. 1.1.4 Distance Functional - for multiple vesicles case

7. 1.2 Outline

8. 1.3 Previous Work and Our Contributions

9. 2 Mathematical Background

10. 2.1 Tangential Calculus

11. 2.2 Shape Differential Calculus

12. 3 Model and Problem Setting

13. 3.1 Willmore energy with Area Constraint

14. 3.1.1 The Model 1 and the Functional $J_{1}(\Gamma)$

15. 3.1.2 Problem Setting of Model 1

16. 3.2 Willmore Energy with Area Constraint and Obstacles

17. 3.2.1 The Model 2 and the Functional $J_{2}(\Gamma)$

18. 3.2.2 Problem Setting of Model 2

19. 3.3 Multiple Vesicles Case

20. 3.3.1 The Distance Functional

21. 3.3.2 The Model 3 and the Functional $J_{3}(\Gamma)$

22. 3.3.3 Problem Setting of Model 3

23. 4 Numerical Schemes

24. 4.1 Time Discretization and Equation Split

25. 4.1.1 Time Discretization

26. 4.1.2 Split

27. 4.2 Finite Elements

28. 4.3 Discrete Problems

29. 4.3.1 Linearization of $dH$ and $dD$ in FEM

30. 4.4 Algorithm

31. 4.4.1 Area Constraint

32. 4.4.2 Full Algorithm

33. 5 Numerical Examples

34. 5.1 Examples

35. 5.2 Golgi Stacks Examples

36. Bibliography

Abbreviations and Notations

Chapter 1 Introduction

The thesis project is a mathematical modeling on the shape evolution of Golgi stack. This project was proposed by Prof. Kang from Life Science Department, and the final goal of this project is to mimic the growing process of Golgi stack so as to give a potential explanation to some special properties of the observed shapes of Golgi stack mathematically. The following image is shown to give a first impression of the Golgi stack [1].

A Golgi stack consists of multiple layers. Each layer is a Golgi cisternae. The cisternae in different stages has different morphological performances. Those stages are called cis-, med- and trans-. For the cis-Golgi, the cisternal assembly is in process, but in the trans-Golgi, the assembly is finished. Our model mainly concerns with trans-Golgi cisternae, for which we can view each layer(cisternae) as a vesicle with closed surface without topological changes.

Our math model is constructed by forming a series of geometric evolution equations. Suppose $\Gamma$ is the surface of a vesicle, which represents one Golgi cisternae. There are four shape functionals considered in our model defined on $\Gamma$. First, the Willmore energy $W(\Gamma)=\int_{\Gamma}h^{2}$, where $h$ is the mean curvature. It is believed that the shapes of the biomembranes are closely relevant to the elastic energy. Because Willmore energy is equivalent to the elastic energy under certain condition, which is explained in Section 1.1.1, we regard the Willmore energy as the dominant energy in our model of Golgi stack. Second, the area functional $A(\Gamma)=\int_{\Gamma}1$. The model include this functional to serve for the surface area constraint. Because the number of molecules is believed to be fixed in our model, the surface area of the shape is conserved. Third, a heaviside functional $H(\Gamma)=\int_{\Gamma}1_{\bf B}$, where $1_{\bf B}$ is an indicator function. To mimic the barriers (the intercisternal elements for Golgi stack) above and belong each Golgi layer, we form this new functional by integrating the heaviside function. Note that the set ${\bf B}$ represents the region of the barriers/obstacles. Fourth, a distance functional $D(\Gamma)$. When multiple vesicles are considered as a whole, the intersection of them should be avoided. Hence, this distance functional is included to handle the relations among those objects. We couple these functionals to form the model. Then we implement it by parametric finite element method, using Matlab.

The followings are the main results of the Golgi simulation. First, for a single layer of trans-cisternae, we mimic the evolution process of the cisternae when the protein vesicles bud from the marginal part of the cisternae. The numerical result explains that it is the barriers placed above and below the cisternae that inhibit the vertical expansion of the cisterae. Second, it is revealed from our result that the swelling of the marginal regions of the trans-cisternae thins its central part to decrease the elastic energy. Third, we also mimic the shape evolution of Golgi with multiple cisternae stacked upon each other. All these simulations are done by mathematical computation, which may give a potential explanation to the mechanism of the Golgi stack.

What’s more, the model can also be extended to more general cases without biological meanings. In Section 5.1, many examples of the applications of this model are demonstrated.

1.1 Overall Introduction and Motivation

In this section, I describe the construction of the math models with explanations on the corresponding biological properties of Golgi stack. From now on, we use the geometric surface $\Gamma\subset\mathbb{R}^{k}$ to describe the surface of a single Golgi cisternae, which is assumed to be a closed vesicle without topological changes. For the case including multiple layers of Golgi cisternaes in one model, we use the family $\{\Gamma_{i}\}_{i=1}^{M}$, where each layer of Golgi cisternae is represented by a surface $\Gamma_{i}$ for some $i$. The followings are four important functionals related to the thesis.

1.1.1 Willmore Energy

First, based on a well-accepted fact suggested in some previous works [2, 3, 4, 5, 6] that the biomembranes are closely related to elastic energy, we consider the elastic energy as the dominated energy to be minimized in our model. We include this energy because $\Gamma$ is also a surface of biomembranes. The elastic energy, also named bending energy, to the lowest order, take the form

[TABLE]

where the integral $\int_{\Gamma}\,dS$ is taken over the surface $\Gamma$, $K_{b}$ is the bending rigidity with respect to mean curvature and $K_{G}$ is the bending rigidity with respect to Gaussian curvature. $h$ and $G$ are the mean and Gaussian curvature respectively defined as $h:=\frac{1}{2}(C_{1}+C_{2}),\ G:=C_{1}C_{2}$. $C_{1}$ and $C_{2}$ are the principle curvatures. Assume that $\Gamma$ is a closed surface. Gauss-Bonnet Theorem [7, Ch. 8] tells,

[TABLE]

where $\chi(\Gamma)$ is the Euler characteristic of $\Gamma$, which is topological invariant. Hence, we only need to consider the Willmore energy [8]

[TABLE]

when $\Gamma$ is a closed surface without topological changes. In summary, to find the optimal shape of Golgi cisternae by minimizing the Willmore energy is the first main point in our model.

1.1.2 Area Constraint

Second, in many cell models, the surface area of a cell is set to be fixed [2, 3]. Agreed by Prof. Kang (Life Science Department), we assume the surface area of each trans-Golgi cisternae $\Gamma$ is conserved. This is based on the fact that the cisternal assembly is completed for tran-cisternae, so the number of molecules of the membrane surface is assumed to be invariant. To enforce this constraint, we consider the functional $A(\Gamma)$ defined by

[TABLE]

which indicates the area of the surface $\Gamma$. Utilizing this functional, we impose the area constraint into the model.

1.1.3 Barrier Functional

Third, we want to include the inter-cisternal elements of Golgi stack, which is biologically relevant to the membrane stacking, into our model. The inter-cisternal elements may limit the expansion of Golgi cisternae in the lateral direction. It serves as the obstacles/barriers placed above and below each Golgi cisternae. To model these restraining factors, we construct a shape functional $H(\Gamma)$ by the integration of an indicator function $1_{\bf B}:\Gamma\rightarrow\{0,1\}$, where ${\bf B}$ is a subset of $\mathbb{R}^{n}$ indicating the region of the obstacles. We define the functional $H(\Gamma)$ as

[TABLE]

where ${\bf x}$ is the identity on $\Gamma$. We want to include this functional to our model so that the surface $\Gamma$ gets hard to cross the region indicated by the set ${\bf B}$. Detailed explanation of the usage of this functional is given in Section 3.2.

1.1.4 Distance Functional - for multiple vesicles case

Fourth, Golgi stack consists of multiple layers of Golgi cisternaes. When multiple vesicles $\{\Gamma_{i}\}_{i=1}^{n}$ are placed in one model, the interaction among the vesicles should be considered. In other words, the vesicles $\Gamma_{i}$ should not cross each other, and even the repulsion between the vesicles should be taken into consideration, because the lipid bilayer structure of the bio-membranes could cause repulsion when the vesicles approach each other. In this case, we consider the distance among the cisternaes. We also use it to construct a functional $D(\Gamma)$. The detailed formulas and explanations are demonstrated in Section 3.3. In summary, $D(\Gamma)$ is built to control the multiple vesicles case.

1.2 Outline

• •

Chapter 2: The preliminary definitions and lemmas related to the thesis are stated. Most of them are in the field of Differential Geometry. Though we did not work on the theories of Geometry, the theoretical results worked by the predecessors are important for us to construct the numerical algorithm.

• •

Chapter 3: The detailed constructions of three models are illustrated. Model 2 and Model 3 are main contributions of this thesis. The three models are introduced in Section 3.1, 3.2 and 3.3 respectively. In each of the section, the motivation of building the model, the idea of the model and the detailed problem setting of the model are stated.

• •

Chapter 4: In the first part of this chapter, the detailed explanation of the time discretization and the space discretization (the construction of the mesh) is given. In the second part, the linearization of some nonlinear functions are illustrated, and then the fully discretized weak formulas for the discrete problems are written for each model introduced in Chapter 3. Finally, the full algorithm is given.

• •

Chapter 5: The first part of this section gives some numerical examples without biological meanings, only to demonstrate the models and to see the conservation of $A(\Gamma)$ and the decrease of the energy. The second part gives some simulations of Golgi stack, including single cisternae case and multi-cisternae case. Some possible explanations to the morphological properties of the Golgi are given, according to the comparison of the numerical results and the observed Golgi.

1.3 Previous Work and Our Contributions

The mathematical study of the shape of biomembranes is introduced by the Helfrich model in 1970s [2, 3], which is a model aiming to study the equilibrium shape of vesicles dominated by elastic energy (or called bending energy). After that, further works on this topic have been done [5, 6, 9, 11, 10, 12, 13, 14, 15], including theoretical analysis and numerical implementation. Many of them applied FEM [11, 10]. Besides, other methods were also studied, for instance, finite difference method [13], level set method [14] and discrete Willmore flow method [15]. Previous works give us inspirations. For example, the Lagrange Multiplier Method is commonly used for the area constraint. Also, previous studies provided excellent formulas for the shape derivative of Willmore energy.

In the thesis, I introduce three models in Section 3.1, 3.2 and 3.3 respectively, and then these models are implemented by FEM with the algorithm stated in Chapter 4. In Section 3.1, it brings out a model dominated by the Willmore energy with conservation of the surface area. This is not a new model developed and solved by us. The works mentioned above [9, 11, 10, 12] have studied this model. Some of them only includes the area constraint [9, 11]. Some also take the volume into account [10, 12]. If one includes both the area and volume constraints, the numerical results could explain the concave shape of the blood cells in a numerical way [2, 10].

Besides, Model 2 (Section 3.2) and Model 3 (Section 3.3) are presented in Chapter 3. These two models are constructed by considering the properties of the Golgi stack and these are the main contribution of our work. Model 2 is an extension of Model 1 by adding some obstacles into the problem. This is inspired by the existence of some biological elements which locate above and below each Golgi cisternae and may confine the vertical extension of the cisternae. Besides, as those confining elements could be moving, we also demonstrate the examples of the moving obstacles. Model 3 is designed for the Golgi stack of multiple cisternaes as a whole. These cisternaes could not cross each other, and even the repulsion of their surface should also be considered. Hence, our method could handle the problem with multiple objects. In summary, Model 2 and 3 are newly formed by us, and we also implement the models by FEM. The numerical results are shown in Chapter 5.

Lastly, by comparing our numerical simulation of the Golgi trans-cisternae to the observed Golgi, some possible explanations are drawn on the morphological performance of the trans-cisternae.

Chapter 2 Mathematical Background

The theoretical background of this research is introduced in this chapter. It contains many theories in the field of Differential Geometry. Without the theoretical results of Geometry worked by the predecessors, the numerical methods will be hard to implement. In this chapter, $\Gamma$ denotes a hyperspace embedded in $\mathbb{R}^{k}$, which is piecewise smooth. The mathematical concepts and theories defined on $\Gamma$, which are relevant to this thesis research, are presented.

2.1 Tangential Calculus

Definition 2.1.1**.**

**(Tangential Gradient)

**The tangential gradient of a function $\phi\in C^{1}(\Gamma;\mathbb{R})$ is defined as

[TABLE]

where $\tilde{\phi}$ is a smooth extension to $\phi$ such that $\tilde{\phi}|_{\Gamma}=\phi$.

Remark 2.1.2*.*

The value of $\nabla_{\Gamma}\phi$ is independent on the extension function $\tilde{\phi}$ chosen.

With the above definition, we can write down the corresponding tangential Jacobian matrix for a vector function $\bm{\phi}\in C^{1}(\Gamma;\mathbb{R}^{d})$:

[TABLE]

Definition 2.1.3**.**

**(Tangential Divergence)

**The tangential divergence of a function $\bm{\phi}\in C^{1}(\Gamma;\mathbb{R}^{d})$ is defined as

[TABLE]

where $\text{tr}(\cdot)$ represents the trace of the matrix.

Definition 2.1.4**.**

**(Tangential Laplace-Beltrami operator)

**The tangential Laplace of a function $\phi\in C^{2}(\Gamma;\mathbb{R})$ is defined as

[TABLE]

2.2 Shape Differential Calculus

Consider a domain $\Gamma\subset{\mathbb{R}^{k}}$ and a functional defined on $\Gamma$ in the following form

[TABLE]

Consider a family of transformations of $\Gamma$, $\{\Gamma_{t}\}_{t\in[0,T]}$. Denote the transformation by $\mathcal{T}_{t}$ such that $\Gamma_{t}=\mathcal{T}_{t}(\Gamma)$ and $\Gamma_{0}=\Gamma$. We assume that $\mathcal{T}_{t}$ is a diffeomorphism from $\Gamma$ to $\Gamma_{t}$ (see [17, Ch. 5] for the definition of diffeomorphism). Denote by $\Gamma_{all}$ the domain containing $\overline{\Gamma_{t}}$ for all $t\in[0,T]$. Let $\bm{v}$ be the velocity field associated with the transformation $\mathcal{T}_{t}$, then we can arrive at the following definition.

Definition 2.2.1**.**

(shape derivative) Suppose $J$ is shape differentiable at $\Gamma$ (see [17, Ch. 5] for the definition). The shape derivative $dJ$ at $\Gamma$ according to the direction $\bm{v}$ is defined as

[TABLE]

where the transformation $\Gamma_{t}=\mathcal{T}_{t}(\Gamma)$ is associated with $\bm{v}$.

Here are two examples of the formulas of shape derivatives of the functional $W(\Gamma)=\int_{\Gamma}h^{2}\,d\Gamma$, where $h$ is the mean curvature. And the detailed derivations of these formulas are provided in [10] and [11]. Note that ${\bf h}:=h\bm{\nu}$.

[TABLE]

[TABLE]

where $D(\phi)_{ij}=[\nabla_{\Gamma}\phi_{j}]_{i}+[\nabla_{\Gamma}\phi_{i}]_{j}$.

Recall the functional

[TABLE]

Now suppose the function value of $\phi$ only depends on ${\bf x}$. In other words, the function value can be represented by $\phi({\bf x})$. According to [18, Ch. 2], the following lemma is applicable.

Lemma 2.2.2**.**

Suppose $\phi\in W_{1}^{2}(\mathbb{R}^{k})$ is independent on the geometry $\Gamma\subset{\mathbb{R}^{k}}$ and $\Gamma$ is of class $C^{2}$. Then in the direction $\bm{v}$,

[TABLE]

For example, this lemma can be applied to the functional $A(\Gamma)=\int_{\Gamma}1\,d\Gamma$. Here $\phi=1$ is clearly independent of $\Gamma$. One can simply apply Lemma 2.1.6 to obtain the formula

[TABLE]

Since $\int_{\Gamma}\text{div}_{\Gamma}\bm{v}\,d\Gamma=\int_{\Gamma}\bm{v}\cdot{\bf h}+\int_{\partial\Gamma}\bm{v}\cdot\bm{\nu}_{s}$ (see [19]), where $\bm{\nu}_{s}$ represents the conormal vector, and the surface $\Gamma$ we consider do not have boundary. Hence, we arrive at the following equation that we widely use in our work:

[TABLE]

Chapter 3 Model and Problem Setting

3.1 Willmore energy with Area Constraint

With the functionals introduced above, we can define a functional $J(\Gamma)$ dependent on the specific problem we are working on, and this functional $J(\Gamma)$ functions as the dominant energy that we want to minimize on the surface $\Gamma$. To track the motion of $\Gamma$ dominated by the energy $J(\Gamma)$, a typical way is to define a geometric evolution equation using the shape derivative $dJ$. Hence, the main idea of our numerical method to solve the shape evolution $\Gamma(t)$ is to find the velocity $\bf v$, which satisfies the equation

[TABLE]

where $(Hil(\Gamma),\langle\cdot,\cdot\rangle)$ is a Hilbert space of functions defined on $\Gamma$.

Consequently, we can define different models by constructing different energies $J(\Gamma)$. In the following sections, I formula three sets of problems corresponding to three models, which can be numerically solved by the discrete schemes described in Chapter 4. In the following three sub-sections, each model is illustrated with detailed equations.

3.1.1 The Model 1 and the Functional $J_{1}(\Gamma)$

The model 1 is a basic model. It considers only one vesicle $\Gamma$. The main idea of Model 1 is to minimize the Willmore energy under the condition that the surface area of $\Gamma$ is fixed. Hence, it does with the functionals $W(\Gamma)=\frac{1}{2}\int_{\Gamma}h^{2}\,d\Gamma$ and $A(\Gamma)=\int_{\Gamma}1\,d\Gamma$. They are briefly discussed in Section 2.1.1 and Section 2.1.2. The biological reason of the conservation of the area is also given in Section 2.1.2.

[TABLE]

where $\Gamma_{0}$ is a given initial shape. The confinement $A(\Gamma)=A(\Gamma_{0})$ is imposed to the problem by using a multiplier $\lambda\in\mathbb{R}$. Then $J_{1}(\Gamma)$ is formulated as

[TABLE]

We aim to find the optimal $\lambda$ and $\Gamma$.

3.1.2 Problem Setting of Model 1

First, we define

[TABLE]

and ${\bf u}:G_{T}\rightarrow\mathbb{R}^{n}$ by

[TABLE]

for all $x\in\Gamma(t)$ and $t\in[0,T]$. Hence, we can regard ${\bf u}(\cdot,t)=\text{Id}_{\Gamma(t)}$. The goal is now to minimize the functional $J_{1}(\lambda,\Gamma)$. Hence, it is a typical way to form an evolution equation in the form of (3.1). Consequently, we have the following problem setting.

Problem 3.1.1**.**

**(Willmore Flow with Area Constraint, Weak Form)

**For a given initial shape $\Gamma_{0}=\Gamma(0)\subset\mathbb{R}^{n}$, find the multiplier $\lambda:[0,T]\rightarrow\mathbb{R}$ and the function ${\bf u}:G_{T}\rightarrow\mathbb{R}^{n}$ according to the family of surface $\{\Gamma(t)\}_{t\in[0,T]}$, such that on the time interval $t\in[0,T]$,

[TABLE]

[TABLE]

for all test function $\phi$. The function space of $\phi$ will be later discussed.

The shape derivatives of the above functionals can be computed using the following equations

[TABLE]

or

[TABLE]

and

[TABLE]

The details of the first two equations can be found in [10] and [11]. The third equation is simpler, so its derivation is given in the Preliminary Section.

Some of the notations used in the above equations are explained as follows:

• •

${\bf h}:=h\bm{\nu}$ is the vector form of mean curvature $h$ on the direction of the outer unit normal vector $\bm{\nu}$.

• •

$\text{Id}={\bf u}(\cdot,t)$ is the identity on $\Gamma(t)$.

• •

$D$ is a symmetric tensor defined by $D(\phi)_{ij}=(\nabla_{\Gamma})_{i}\phi^{j}+(\nabla_{\Gamma})_{j}\phi^{i}$.

The equation (3.7) and (3.9) are implemented for the Model 1.

This Model 1 is not a model firstly produced and solved by us. Instead, many works [9, 11, 10, 12] have been done on the study of this model, of which many also consider the volume constraint [10, 12]. In summary, Model 1 is a fundamental model, which has been studied for many years. However, it catches an important property of the bio-membranes, that the membranes is closely relevant to the elastic energy. Also, it imposes the idea that the number of molecules of the membranes is fixed such that the surface area of the membrane is conserved. If one includes one more constraint, the volume constraint, to Model 1, the numerical results can explain the concave shape of the blood cells mathematically [2, 10]. However, since our work is triggered by the shape evolution of Golgi stack, instead of considering the volume constraints, we consider other special properties of the Golgi stack and build our own models. The coming two sections present two models of Golgi stacks. Model 2 is designed for a singer layer of Golgi cisternae and Model 3 is designed for the Golgi stack of multiple cisternaes as a whole.

3.2 Willmore Energy with Area Constraint and Obstacles

3.2.1 The Model 2 and the Functional $J_{2}(\Gamma)$

Model 2 is an extension to Model 1 by considering the existence of some obstacles/barriers in the problem. This is motivated by a property of Golgi cisternae, which may be confined by the intercisternal elements above and below each Golgi cisternae layer. (See Section 2.1.3). In Model 2, the obstacles are considered and estimated by the functional

[TABLE]

defined in Section 2.1.3, where ${\bf B}$ represents the region of the obstacles/barriers. One can easily observe that $H(\Gamma)$ will take nonzero values only if $\Gamma\cap{\bf B}\neq\emptyset$. Hence, if we add this functional $H(\Gamma)$ to the energy that we want to minimize, it will be hard for the shape $\Gamma$ to touch and cross the region ${\bf B}$ so as to avoid the increase of the total amount of the energy.

To conclude, in Model 2, we want to do the same optimization as in Model 1, but also to include some obstacles indicated by the set ${\bf B}\subset\mathbb{R}^{n}$. Normally, this model can be explained by the following problem:

[TABLE]

The resulted functional $J_{2}$ can then be formulated as

[TABLE]

which is the augmented energy to be minimized in Model 2.

Remark 3.2.1*.*

The constant $\alpha\in\mathbb{R}$ is a weight of the functional $H(\Gamma)$ to control the impact of $H$. The dominated energy of our model should be the Willmore energy $W(\Gamma)$. $H(\Gamma)$ is only a constraint. Hence, we don’t want the functional $H(\Gamma)$ to dominate the whole energy $J_{2}(\Gamma)$.

3.2.2 Problem Setting of Model 2

Use the same notation as those in Section 3.1.2, we can form the following weak problem to minimize $J_{2}(\lambda,\Gamma)$.

Problem 3.2.1**.**

**(Willmore Flow with Area Constraint and Obstacles, Weak Form)

**Suppose that $\alpha\in\mathbb{R}$ is fixed as a weight coefficient. Now given $\Gamma_{0}=\Gamma(0)\subset\mathbb{R}^{n}$, find $\lambda:[0,T]\rightarrow\mathbb{R}$ and the function ${\bf u}:G_{T}\rightarrow\mathbb{R}^{n}$ such that on the time interval $t\in[0,T]$,

[TABLE]

[TABLE]

for all test function $\phi$.

The calculation of $dH$ is discussed as follows.

Since the indicator function $1_{\bf B}:\Gamma\rightarrow\{0,1\}$ is discontinuous, the implementation of it using FEM is not applicable. In the finite element method, we choose a smooth version of the indicator. More precisely, assume that ${\bf B}$ represents a very regular shape. Then, $1_{\bf B}$ can be written as a composition of the heaviside function

[TABLE]

Remark 3.2.2*.*

For example, if ${\bf B}=[a,b]\times([m,\infty]\cup[-\infty,-m])\subset\mathbb{R}^{2}$, then it can be written as a composition of $heav(x)$ by

[TABLE]

Though $heav(x)$ is still a discontinuous step function, we use the smooth approximation

[TABLE]

of it. Now a smooth version of $1_{\bf B}$ is obtained.

With the smoothness of $1_{\bf B}$, we can now derive the formula for $dH$. By lemma 2.2.2, the formula of $dH(\Gamma)$ is obtained:

[TABLE]

Remark 3.2.3*.*

Usually, in the math models, ${\bf B}$ is some fixed obstacles and independent of time. However, inspired by the hypothesis mentioned by Prof. Kang, that the intercisternal elements (such as Golgi matrix, which works as the constraining factor in our math model on Golgi stacks) maintain the same distance with the membrane since some of its components are embedded in the membrane. Hence, it could be more realistic to keep the distance between the membrane and the barriers when we mimic the growing process of Golgi cisternae. Nevertheless, since the shape of the cisternae $\Gamma(t)$ evolves with time, the position of the barriers also need changesx to keep the distance. Based on this, we also make some numerical experiments for the moving obstacles ${\bf B}(t)$. These examples can be found in Chapter 5.

3.3 Multiple Vesicles Case

3.3.1 The Distance Functional

As succinctly introduced in Section 2.1.4, the Golgi stack consist of multiple layers of Golgi cisternaes, denoted by $\Gamma_{1},\Gamma_{2},\cdots,\Gamma_{n}$ in our math model. The family $\{\Gamma_{i}\}_{i=1}^{n}$ is aimed to be modeled on the whole, at the same time the vesicles $\Gamma_{i}$ should not cross or even should repulse from each other. We applied the following function with Euclidean norm $\|\cdot\|$ to measure the distance between the vesicles:

[TABLE]

where $V_{i}$ denotes the $i$-th vesicle. Equivalently, $d(x)$ can be defined as

[TABLE]

Using this measurement of distance, we define the following functional

[TABLE]

where $x$ is the identity on $\Gamma$. Note that $\Gamma$ only represents the surface of one single vesicle $V_{i}$.

3.3.2 The Model 3 and the Functional $J_{3}(\Gamma)$

Well-prepared with the above functional $D(\Gamma)$, we can now extend the single vesicle case - Model 2 to the multiple case - Model 3. Let $\Gamma=\partial V_{i}$ for some $i$. For each $i$, we do the following problem:

[TABLE]

where $\alpha,\beta\in\mathbb{R}$ are the weight coefficients for $H(\Gamma)$ and $D(\Gamma)$ respectively. To translate this optimization problem: in fact, we are minimizing the Willmore energy $W(\Gamma)$, under the condition that, firstly, $\Gamma$ is hard to tough the region $B$ with a weight $\alpha$ and secondly, $\Gamma$ is hard to get very closed to other vesicles with a weight $\beta$.

The associated functional $J_{3}$ corresponding to this problem is given by

[TABLE]

which is the augmented energy to be minimized in Model 3.

Remark 3.3.1*.*

The Model 3 is also newly formed by us. It is motivated by the component of the multiple layers of Golgi. It could be extended to other problems which include multiple objects.

3.3.3 Problem Setting of Model 3

Similar to what we do for the above two models, we can now form the following weak problem to minimize the functional $J_{3}(\lambda,\Gamma)$.

Problem 3.3.1**.**

**(Willmore Flow with Area Constraint and Obstacles - applied to Multiple Vesicles Case, Weak Form)

**Suppose that $\alpha,\beta\in\mathbb{R}$ are fixed as a weight coefficient for the functional $H$ and $D$ respectively. Given an initial shape $\Gamma_{0}=\Gamma(0)\subset\mathbb{R}^{n}$, find the multiplier $\lambda:[0,T]\rightarrow\mathbb{R}$ and the function ${\bf u}:G_{T}\rightarrow\mathbb{R}^{n}$ according to the family of surfaces $\{\Gamma(t)\}_{t\in[0,T]}$, such that on the time interval $t\in[0,T]$,

[TABLE]

[TABLE]

for all test function $\phi$.

The formulas of the shape derivatives of $W(\Gamma)$, $H(\Gamma)$ and $A(\Gamma)$ are clearly explained in the previous sections, so I only demonstrate the calculation of $dD(\Gamma;\phi)$ in this section. Since our model is implemented by FEM at the end, I estimate the shape derivative of $D(\Gamma)$ in a tricky way. First, for any point $x^{*}\in\Gamma_{i}$, consider an open ball $S=B(x^{*};r)$ centered at $x^{*}$ with radius $r$ and define a function $d^{*}:S\rightarrow\mathbb{R}^{+}$ by

[TABLE]

Since $1/d^{*}$ is independent on the geometry, we can apply Lemma 2.2.2 and then obtain the shape derivative of $D^{*}(S\cap\Gamma):=\int_{s\cap\Gamma}f^{*}(x)\,d\Gamma$.

[TABLE]

Applying this formula, then we can approximate the shape derivative of $D(\Gamma)$ by piecewise implementation in FEM.

Chapter 4 Numerical Schemes

4.1 Time Discretization and Equation Split

4.1.1 Time Discretization

The model is implemented on the time-interval $[0,T]$, though the final time $T$ can be chosen dependent on the stopping criteria set in the algorithm. Let

[TABLE]

where $t_{0}=0,\ t_{N}=T$, be a partition of the interval $[0,T]$.

Remark 4.1.1*.*

About Time Adaptivity: If one wants to obtain a more effective algorithm to reach the optimization shape of $\Gamma$ faster, it is more reasonable to make the time step

[TABLE]

adaptive to the mesh size, because this FEM is using a moving-mesh. For me, I simply choose a comparatively small time step $\tau=\tau_{n}$, which is fixed, for convenience. However, there is actually a goodness of using a small time step $\tau$, because we apply the linear approximations (see Section 4.3.1) on some shape functionals in our method and a small time step is beneficial to the linear approximations.

Now the solution we want to find is the family $\{\Gamma(t_{i})\}_{i=1}^{N}$ when $\Gamma(t_{0})$ is given. Recall the function ${\bf u}:G_{T}\rightarrow\mathbb{R}^{n}$ defined by equation (3.4). From now on, we denote the numerical solution by $\Gamma^{n}$ at each time $t_{n}$, which is viewed as the image of ${\bf U}^{n}(\cdot)$

[TABLE]

with given approximation $\Gamma_{0}=\Gamma^{0}\approx\Gamma(t_{0})$.

4.1.2 Split

From equation (3.5) and (3.7), one can see that, if we solve ${\bf u}$ directly, the order of the differential equation is high. To reduce this order, many previous work chose to split the formula [10]. In our work, we solve a pair of unknown $({\bf V}^{n+1},{\bf H}^{n+1})$ first, and then update ${\bf U}^{n+1}$. ${\bf V}^{n+1}$ is defined as

[TABLE]

By applying the equation $-\Delta_{\Gamma}x={\bf h}$ [19, Page 390] and the above equation (4.2), the following identity for ${\bf H}^{n+1}$ is obtained

[TABLE]

and hence the following weak formula

[TABLE]

for all test function $\phi$. The function space of the test function will be discussed in the coming section.

In summary, with equation (4.4), we can solve the problems stated in Chapter 3 by solving the pair of unknowns $({\bf V}^{n+1},{\bf H}^{n+1})$ first and then update ${\bf U}^{n+1}$ by equation (4.2).

4.2 Finite Elements

The following are some notations and definitions used in this section.

• •

$\mathbb{R}^{k}$ is the space containing $\Gamma$.

• •

$\mathbb{R}^{k-1}$ is the parametrization space.

Polyhedral Approximation $(\overline{\Gamma}_{h},\overline{\mathcal{T}}_{h},\overline{K},V_{\overline{K}})$

• •

$\overline{\Gamma}_{h}=\bigcup_{\overline{K}\in\overline{\mathcal{T}}_{h}}\overline{K}$ is a polyhedral approximation of $\Gamma$, where $\overline{\mathcal{T}}_{h}$ is the triangulation of $\overline{\Gamma}_{h}$ and the vertices of $\overline{\Gamma}_{h}$ lie on $\Gamma$.

• •

$\overline{K}$ is a $(k-1)$-simplex in $\mathbb{R}^{k}$ with its $k$ vertices $\{v_{i}\}_{i=1}^{k}$.

• •

$V_{\overline{K}}=\{v=\sum_{i=1}^{k}c_{i}v_{i};\ \ \sum_{i=1}^{k}c_{i}=1,c_{i}\in\{0,\frac{1}{2},1\}\}$ is the vertex set that we use in FEM, which includes the vertices of $\overline{K}$ and the mid-points of each edge of $\overline{K}$.

Polynomial Approximation $(\Gamma_{h},K,V_{K})$

• •

$\Gamma_{h}$ is the image of a function $p$ defined on $\overline{\Gamma}_{h}$ such that $p|_{\overline{K}}$ is a polynomial of degree $\leq 2$.

• •

Denote $K=p(\overline{K})$. Then we have $\Gamma_{h}=\bigcup K$.

• •

$V_{K}=\{p(v);\ v\in V_{\overline{K}}\}$.

• The function values of $p(v)$ indicate the position of the vertices $V_{K}$. The values of $p(v)$ is obtained by following the rules:

If $v\in\cup_{i=1}^{k}\{v_{i}\}$ (i.e. $v$ is a vertex of the $(k-1)$-simplex $\overline{K}$), then $p(v)=v$. Since $v_{i}$ are on $\Gamma$, then $p(v_{i})$ are also on $\Gamma$. 2. 2.

If $v\in V_{k}\setminus\cup_{i=1}^{k}\{v_{i}\}$ (i,e, $v$ is a mid-point of an edge of $\overline{K}$), then $p(v)$ is an orthogonal projection onto $\Gamma$. Hence, $p(v)$ also lies on $\Gamma$.

To make sure that the set $V_{K}^{n+1}$ is the orthogonal projection of $V_{\overline{K}}$ on to $\Gamma_{h}^{n+1}$, we adjust the position of the vertices whenever we obtain a new set $V_{K}^{n+1}$ from $V_{K}^{n}$ by solving the numerical problem.

By the above construction, one can conclude that all the points in $V_{K}$ are on $\Gamma$. Besides, the existence and uniqueness of the polynomial $p|_{K}$ satisfying the above two rules are proved by [16, Theorem 2.2.1].

Reference Element

• •

$K_{ref}$, a $(k-1)$-simplex in $\mathbb{R}^{k-1}$. We define the standard reference $K_{ref}$ as the convex envelope with vertices $\cup_{i=1}^{k-1}\{{\bf e}_{i}\}\cup\{{\bf 0}\}$.

• –

In $\mathbb{R}^{1}$, $K_{ref}=[0,1]$.

• –

In $\mathbb{R}^{2}$, $K_{ref}=\{a_{0}(0,0)+a_{1}(0,1)+a_{2}(1,0):\ \forall a_{0},a_{1},a_{2}\geq 0\ \text{s.t.}\ a_{0}+a_{1}+a_{2}=1\}$.

• •

$V_{ref}=\{v=\sum_{i=1}^{k-1}c_{i}{\bf e}_{i};\ \ \sum_{i=1}^{k}c_{i}=1,c_{i}\in\{0,\frac{1}{2},1\}\}$ represents the vertex set of $K_{ref}$.

• For each $(k-1)$-simplex $\overline{K}$ in $\mathbb{R}^{k}$ mentioned above, there exists a bijective mapping $f_{\overline{K}}:K_{ref}\rightarrow\overline{K}$ such that $f_{\overline{K}}$ maps $V_{ref}$ to $V_{\overline{K}}$.

Finite Element Space

• The finite element space defined over the set $\Gamma_{h}$ is defined as

[TABLE]

• At each time step $t_{n+1}$, $\Gamma_{h}^{n}$ is known. Hence, the finite element function space that we consider at the time step $t_{n+1}$ to find $\Gamma_{h}^{n+1}$ is $\mathbb{F}(\Gamma_{h}^{n})$.

4.3 Discrete Problems

With the time and space discretization discussed in the previous sections, we can now rewrite our problems in the discrete forms, which can be implemented.

In Chapter 3, we formulate the three problems, Problem 3.1.1, Problem 3.2.1 and Problem 3.3.1. In each of them, we aim to solve ${\bf u}(\cdot,t)=\text{Id}_{\Gamma(t)}$ for $t\in[0,T]$. After time discretization, we aim to solve $\Gamma^{n}:=\{{\bf U}^{n}(x);\ \forall x\in\Gamma^{n-1}\}$ (equation (4.1)) at each time $t_{n}$. Instead of solving ${\bf U}^{n+1}$ directly from ${\bf U}^{n}$, as discussed in Section 4.1.2, we split the formula by adding the equation (4.4), so we can solve $({\bf V}^{n+1},{\bf H}^{n+1})$ first and then update ${\bf U}^{n+1}$ by equation (4.2). The discretization (Section 4.2) gives us the space $\mathbb{F}(\Gamma_{h}^{n})$ to solve $({\bf V}^{n+1},{\bf H}^{n+1})$. Based on these ideas, we can formulate the problems in Chapter 3 into the following discrete forms.

Problem 4.3.1**.**

**(Discrete Form of Problem 3.1.1.)

**Suppose that $\alpha,\beta\in\mathbb{R}$ are fixed as a weight coefficient for the functional $H$ and $D$ respectively. Given an initial shape $\Gamma_{h}^{0}\approx\Gamma(0)\subset\mathbb{R}^{k}$, find ${\bf V}^{n+1},{\bf H}^{n+1}\in\mathbb{F}(\Gamma_{h}^{n})$ and $\lambda^{n+1}\in\mathbb{R}$ such that for each $n$, $\forall\Phi\in\mathbb{F}(\Gamma_{h}^{n})$,

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

and

[TABLE]

At each time step, $\Gamma_{h}^{n+1}$ is updated by

[TABLE]

Problem 4.3.2**.**

**(Discrete Form of Problem 3.3.1.)

**Given an initial shape $\Gamma_{h}^{0}\approx\Gamma(0)\subset\mathbb{R}^{k}$, find ${\bf V}^{n+1},{\bf H}^{n+1}\in\mathbb{F}(\Gamma_{h}^{n})$ and $\lambda^{n+1}\in\mathbb{R}$ such that for each $n$, $\forall\Phi\in\mathbb{F}(\Gamma_{h}^{n})$,

[TABLE]

[TABLE]

[TABLE]

where the formulations of $dH^{n+1}$ and $dD^{n+1}$ are linearized in the next section. They are specially treated and formulated in an implicit form.

Note that the discrete form of Problem 3.2.1. is the same as Problem 4.3.2. by simply taking $\beta$ to be zero.

4.3.1 Linearization of $dH$ and $dD$ in FEM

Linearization of $dH$

Recall the formula for $dH$ (3.14):

[TABLE]

The function $1_{\bf B}$ with smoothness (by applying exponential functions) is nonlinear. So is $\nabla 1_{\bf B}$. In FEM, we use their linear approximation. They are simply made in the standard way:

[TABLE]

[TABLE]

Here $\mathcal{D}$ is the Fréchet derivative operator. It is well-known that the approximation is good for ${\bf x}$ if it is close enough to ${\bf x_{0}}$, so it is very natural that when ${\bf x}\in\Gamma^{n+1}$ we take the points ${\bf x_{0}}$ from $\Gamma^{n}$. When $\tau$ is small, the approximation could be good. Hence, we now replace ${\bf x_{0}}$ by $x\in\Gamma_{h}^{n}$ and ${\bf x}$ by ${\bf U}_{h}^{n+1}(x)\in\Gamma_{h}^{n+1}$. More specifically, we write the approximation as: $\forall x\in\Gamma_{h}^{n}$,

[TABLE]

[TABLE]

Clearly, in the discretized form, we can replace $\left({\bf U}_{h}^{n+1}(x)-{x}\right)$ by $\left(\tau{\bf V}^{n+1}\right)$ and then obtain the following semi-implicit formula for $dH^{n+1}(\Gamma_{h}^{n};\Phi)$:

[TABLE]

Linearization of $dD$

Recall the formula (3.22) for $dD$ on the neighborhood $S\cap\Gamma$ of a point $x^{*}\in\Gamma$:

[TABLE]

Similarly, we first linearize the functions $\dfrac{1}{d^{*}(x)^{2}}(x-y^{*})$ and $\dfrac{1}{d^{*}(x)}$ in the standard way:

[TABLE]

[TABLE]

Here $\mathcal{H}\left(\dfrac{1}{d^{*}({\bf x_{0}})}\right)$ is the Hessian matrix of the function $\dfrac{1}{d^{*}}$ at the point ${\bf x_{0}}$. Similarly, we take $x_{0}$ to be the points ${\bf x}$ on $\Gamma_{h}^{n}$ and then obtain the following equations for all ${\bf x}\in\Gamma_{h}^{n}$:

[TABLE]

[TABLE]

With this approximation, we construct the formula for $dD^{n+1}(\Gamma_{h}^{n};\Phi)$:

[TABLE]

Note that $d_{K}({\bf x})$ is exactly the distance function $d^{*}({\bf x})$ by taking the point ${\bf y^{*}}$ dependent on each element $K$, denoted by ${\bf y}_{k}$. Understandably, ${\bf y}_{k}$ is taken in the multiple vesicles case where ${\bf y}_{k}$ is a point on one of the vesicles other than $\Gamma_{h}^{n}$ such that it is closest to the element $K$.

4.4 Algorithm

The problem is fully discretized as discussed in the last chapter. Since the Problem 4.3.2. is the full problem when the other two problems can be obtained by taking either $\alpha$ or $\beta$ to be zero, we illustrate the full algorithm for this problem in this chapter. We now have the Problem 4.3.2. discretized and linearized, which make it standard to be solved by FEM, except that the area constraint $A(\Gamma_{h}^{n+1})=A(\Gamma_{h}^{0})$ need to be reconsidered. Hence, we start our discussion from the part of ’area constraint’.

4.4.1 Area Constraint

The method to treat the area constraint is based on the method presented in [10] to compute the Lagrange multiplier $\lambda^{n+1}$, but with some nontrivial difference.

First, recall the equations for Problem 4.3.2. here for convenience:

[TABLE]

[TABLE]

[TABLE]

In their method, they make the shape derivative $dA$ in the above equation to be explicit [10]. That is to say, it becomes $dA^{n}(\Gamma_{h}^{n};\Phi)=\int{\bf H}^{n}\cdot\Phi$ instead of $dA^{n+1}(\Gamma_{h}^{n};\Phi)=\int{\bf H}^{n+1}\cdot\Phi$ as stated in our method. However, when I try to make it explicit as they said and implement the method in Matlab, the program always breaks down. Even if I take $\alpha$ and $\beta$ to be zero, which makes our model very similar to their problem, the numerical results still break down. However, when I use the implicit formula $dA^{n+1}$ as in equation (4.21), the numerical results reveal to be stable and convergent. Hence, we apply the ideas of computing Lagrange Multiplier as stated in [10], but use $dA^{n+1}$. The idea for solving our problem exactly stated as in equations (4.21 - 4.4.2) is explained below.

Rewrite the pair of unknown as $({\bf V}^{n+1},{\bf H}^{n+1})=({\bf V}^{n+1}_{1},{\bf H}^{n+1}_{1})+\lambda^{n+1}({\bf V}^{n+1}_{2},{\bf H}^{n+1}_{2})$. We solve $({\bf V}^{n+1}_{1},{\bf H}^{n+1}_{1})$ and $({\bf V}^{n+1}_{2},{\bf H}^{n+1}_{2})$ separately corresponding to the Problem 4.4.1 and Problem 4.4.2, and then find the Lagrange Multiplier $\lambda^{n+1}$ such that the resulted $\Gamma_{h}^{n+1}$ updated by ${\bf V}^{n+1}={\bf V}^{n+1}_{1}+\lambda^{n+1}{\bf V}^{n+1}_{2}$ satisfies the area constraint.

Problem 4.4.1**.**

Find $({\bf V}^{n+1}_{1},{\bf H}^{n+1}_{1})$ such that $\forall\Phi\in\mathbb{F}(\Gamma_{h}^{n})$,

[TABLE]

[TABLE]

Problem 4.4.2**.**

Find $({\bf V}^{n+1}_{2},{\bf H}^{n+1}_{2})$ such that $\forall\Phi\in\mathbb{F}(\Gamma_{h}^{n})$,

[TABLE]

[TABLE]

The above two problems are standardly formulated to be solved by FEM. Now the values of $({\bf V}^{n+1}_{1},{\bf H}^{n+1}_{1})$ and $({\bf V}^{n+1}_{2},{\bf H}^{n+1}_{2})$ are known and we only need to find the suitable $\lambda^{n+1}$. More specifically, by equation , we solve the Lagrange Multiplier $\lambda^{n+1}\in\mathbb{R}$ by Newton’s Method as a root of the function

[TABLE]

where

[TABLE]

The detailed information about the derivation of the differential of $f_{N}$ and the initial guess for this Newton’s method is given by [10]. The iterative equation is

[TABLE]

where $\mathcal{D}f_{N}(\lambda_{i-1})=\tau\int_{\Gamma^{n}(\lambda)}\text{dir}_{\Gamma^{n}(\lambda)}{\bf V}_{2}^{n+1}$. The initial guess is

[TABLE]

4.4.2 Full Algorithm

Chapter 5 Numerical Examples

In this chapter, we demonstrate various numerical examples to show the various applications of the models (from Model 1 to Model 3). In Section 5.1, the experiments are not aimed to mimic something. We just try different initial shapes to make as many interesting experiments as possible. Also, we plot some graphs to see the decrease of energy. In Section 5.2, the experiments are aimed to mimic the shape properties of the Golgi stacks, so the initial shapes are set goal-oriented. In this section, we provide additional information about the biological motivation of forming those special experiments. Also, we combine those numerical results with the observed image of Golgi cisternae and draw some conclusions based on that.

5.1 Examples

In this section, we present various examples, from the most basic one. The first example is a minimization of the functional $A(\Gamma)$. It is simply the length of $\Gamma$ in $\mathbb{R}^{2}$. This model is simple but contained in every model. Figure 5.1 shows that the ellipse shrinks to a point.

Example 1

Example 2 － Model 1

The second example is the minimization of Willmore Energy $W(\Gamma)$ under the condition that the area (length in 2D) $A(\Gamma)$ is fixed.

Example 3 - Model 1

This is similar to example 2 but it has a more interesting initial shape - the C shape.

Example 4 - Model 2

Now we begin to demonstrate some examples of Model 2, where the barriers/obstacles ${\bf B}$ is included. In all these examples, the dash lines or the grey block areas represent the region of the obstacles ${\bf B}$.

Example 5 - Model 2

This is an interesting example. When the barriers are put on one side of an ellipse, the ellipse will get rid of the bounds and then become a circle. This experiment is pretty realistic.

Example 6 - Model 2

The examples continue in next section. In Section 5.2, more examples of Model 2 and Model 3 are explained. In particular, those examples have their practical meanings applying to the shape evolution of Golgi stacks.

5.2 Golgi Stacks Examples

Detailed Biological Explanations

There are different stages (cis-,med-,trans-) for the Golgi cisternaes. The morphological changes we aim to mimic are those for trans-cisternae in plant cells. Biologists have done some works on the ET analyses of Golgi, which observed that trans-cisternae are less thicker than interior cisternae [20, 21]. Prof. Kang (Life Science Department, CUHK) also provided an image of trans-cisternae in the border cell (Figure 5.7), from which one can draw similar conclusion with the previous biological work.

To explain the observed morphological properties, we model the evolution of the Golgi membrane $\Gamma$ in $\mathbb{R}^{2}$ for simplicity. The idea of our modeling is roughly stated in Chapter 2. More specifically, suppose that the shapes of the cisternal membrane are determined to minimize the elastic energy of the membrane with two barriers lying above and below the cisternae. These barriers are used to model some inter-cisternal elements, so for a single cisternae they just lie above and below. This is how we place the obstacles ${\bf B}$. Besides, based on the condition that cisternal assembly is completed in the med- stage of Golgi, the surface area of the membrane only grows in the med- stage, but for trans-cisternae the number of molecules which consist the membrane is fixed. Hence, the constraint $A(\Gamma)=A(\Gamma_{0})$ applies. Lastly, the maturation process of Golgi in the plant border cells involves lots of synthesis of the biological substance. Naturally, the luminal volume of trans-cisternae should increase, though biologists did not measure it specifically. We do not consider a condition on the volume in our math model, because there is too much uncertainty about the measurement of the increase. However, when the elastic energy decrease, the volume usually increase naturally.

Numerical Examples for Single Cisternae

With these ideas, we use Model 2 to mimic the evolution of a single trans-cisternae with initial shape $\Gamma_{0}$ shown in Figure 5.7(a). The first set of numerical results are demonstrated in Figure 5.7.

Observation and Conclusion

By observing the numerical experiment shown in Figure 5.7, we arrive at the following conclusions.

First, When the obstacles ${\bf B}$ were placed to restrain the vertical expansion of the cisternal membrane $\Gamma$, the growth was limited to the peripheral regions of the cisternae. As a comparison, if we remove the obstacles, the cisternae grows to a peanut shape (as shown in Figure 5.11(b)-(c)) and then it grows isotropically. The cisternae without barriers finally become a sphere (Figure 5.11(h)), which is the optimized shape when the length of $\Gamma$ in $\mathbb{R}^{2}$ is fixed. In conclusion, Figure 5.7 shows an mathematical example that matches the phenomenon that the shapes of trans-cisternae evolve out to the cisternal margin by confining the extension of the central domain and decreasing the elastic energy. Note that naturally it is more stable if the energy is low.

Second, it is intriguing to find that the central domain of the cisternae gets thinner when the elastic energy decreases. It is not an essential phenomenon when the marginal domain is swelling. Figure 5.8 demonstrates a counter example. In this example, the central domain is also limited by the barriers. The marginal part also swells. The volume also increases. The length is also fixed. However, the central domain dose not become thinner. Hence, we here relate the thinning of the central part of the trans-cisternae, which is observed by previous biological works [20, 21], to the decrease of the elastic energy of the cisternal membrane.

A Single Cisternae with Moving Barriers $\bf B$

Usually, in the math models, the position of the obstacles ${\bf B}$ is fixed. However, inspired by the hypothesis mentioned by Prof. Kang, that the intercisternal elements (such as Golgi matrix) maintain the same distance with the membrane since some of its components are embedded in the membrane. Hence, it could be more realistic to keep the distance between the membrane and the barriers.

A Single Cisternae without Barriers

Multiple Layers

The Golgi stacks consist of multiple layers of cisternaes (see Figure 5.12). This fact inspires us to form a model of multiple vesicles $\{\Gamma_{i}\}_{i=1}^{M}$.

The details of the multiple vesicles case are illustrated in Section 3.3. The following example is an implementation of Model 3, with the initial shapes given as a set of parallelly ellipse-like shapes. In Figure 5.13, three cisternaes, named $\Gamma_{1},\ \Gamma_{2}$ and $\Gamma_{3}$ from up to down are considered.

We do not draw conclusions about the Golgi stacks based on the modeling of multiple cisternaes, because the components of the Golgi stack is much more complicated. Not only the different biological properties of different stages of cisternaes make it complicated, but also the mechanism between them. Hence, we do these examples as a try only, though we can still find some similarity between our numerical results and the observed Golgi. For example, in Figure 5.12, the marginal parts of the upper half cisternaes go up, while those of the cisternaes below go down. Our numerical results also reveal this tendency.

The simulation of more layers (5 cisternaes in Figure 5.14) gives similar results as those in Figure 5.13.