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 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 , where each layer of Golgi cisternae is represented by a surface for some . 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 is also a surface of biomembranes. The elastic energy, also named bending energy, to the lowest order, take the form
where the integral is taken over the surface , is the bending rigidity with respect to mean curvature and is the bending rigidity with respect to Gaussian curvature. and are the mean and Gaussian curvature respectively defined as . and are the principle curvatures. Assume that is a closed surface. Gauss-Bonnet Theorem [7, Ch. 8] tells,
where is the Euler characteristic of , which is topological invariant. Hence, we only need to consider the Willmore energy 
when 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 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 defined by
which indicates the area of the surface . 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 by the integration of an indicator function , where is a subset of indicating the region of the obstacles. We define the functional as
where is the identity on . We want to include this functional to our model so that the surface gets hard to cross the region indicated by the set . 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 are placed in one model, the interaction among the vesicles should be considered. In other words, the vesicles 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 . The detailed formulas and explanations are demonstrated in Section 3.3. In summary, is built to control the multiple vesicles case.
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 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 , level set method  and discrete Willmore flow method . 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.
2.1 Tangential Calculus
The tangential gradient of a function is defined as
where is a smooth extension to such that .
The value of is independent on the extension function chosen.
With the above definition, we can write down the corresponding tangential Jacobian matrix for a vector function :
The tangential divergence of a function is defined as
where represents the trace of the matrix.
(Tangential Laplace-Beltrami operator)
The tangential Laplace of a function is defined as
2.2 Shape Differential Calculus
Consider a domain and a functional defined on in the following form
Consider a family of transformations of , . Denote the transformation by such that and . We assume that is a diffeomorphism from to (see [17, Ch. 5] for the definition of diffeomorphism). Denote by the domain containing for all . Let be the velocity field associated with the transformation , then we can arrive at the following definition.
(shape derivative) Suppose is shape differentiable at (see [17, Ch. 5] for the definition). The shape derivative at according to the direction is defined as
where the transformation is associated with .
Recall the functional
Now suppose the function value of only depends on . In other words, the function value can be represented by . According to [18, Ch. 2], the following lemma is applicable.
Suppose is independent on the geometry and is of class . Then in the direction ,
For example, this lemma can be applied to the functional . Here is clearly independent of . One can simply apply Lemma 2.1.6 to obtain the formula
Since (see ), where represents the conormal vector, and the surface we consider do not have boundary. Hence, we arrive at the following equation that we widely use in our work:
3.1 Willmore energy with Area Constraint
With the functionals introduced above, we can define a functional dependent on the specific problem we are working on, and this functional functions as the dominant energy that we want to minimize on the surface . To track the motion of dominated by the energy , a typical way is to define a geometric evolution equation using the shape derivative . Hence, the main idea of our numerical method to solve the shape evolution is to find the velocity , which satisfies the equation
where is a Hilbert space of functions defined on .
Consequently, we can define different models by constructing different energies . 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
The model 1 is a basic model. It considers only one vesicle . The main idea of Model 1 is to minimize the Willmore energy under the condition that the surface area of is fixed. Hence, it does with the functionals and . 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.
where is a given initial shape. The confinement is imposed to the problem by using a multiplier . Then is formulated as
We aim to find the optimal and .
3.1.2 Problem Setting of Model 1
First, we define
for all and . Hence, we can regard . The goal is now to minimize the functional . Hence, it is a typical way to form an evolution equation in the form of (3.1). Consequently, we have the following problem setting.
(Willmore Flow with Area Constraint, Weak Form)
For a given initial shape , find the multiplier and the function according to the family of surface , such that on the time interval ,
for all test function . The function space of will be later discussed.
The shape derivatives of the above functionals can be computed using the following equations
Some of the notations used in the above equations are explained as follows:
is the vector form of mean curvature on the direction of the outer unit normal vector .
is the identity on .
is a symmetric tensor defined by.
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
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
defined in Section 2.1.3, where represents the region of the obstacles/barriers. One can easily observe that will take nonzero values only if . Hence, if we add this functional to the energy that we want to minimize, it will be hard for the shape to touch and cross the region 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 . Normally, this model can be explained by the following problem:
The resulted functional can then be formulated as
which is the augmented energy to be minimized in Model 2.
The constant is a weight of the functional to control the impact of . The dominated energy of our model should be the Willmore energy . is only a constraint. Hence, we don’t want the functional to dominate the whole energy .
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 .
(Willmore Flow with Area Constraint and Obstacles, Weak Form)
Suppose that is fixed as a weight coefficient. Now given , find and the function such that on the time interval ,
for all test function .
The calculation of is discussed as follows.
Since the indicator function 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 represents a very regular shape. Then, can be written as a composition of the heaviside function
For example, if , then it can be written as a composition of by
Though is still a discontinuous step function, we use the smooth approximation
of it. Now a smooth version of is obtained.
With the smoothness of , we can now derive the formula for . By lemma 2.2.2, the formula of is obtained:
Usually, in the math models, 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 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 . 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 in our math model. The family is aimed to be modeled on the whole, at the same time the vesicles should not cross or even should repulse from each other. We applied the following function with Euclidean norm to measure the distance between the vesicles:
where denotes the -th vesicle. Equivalently, can be defined as
Using this measurement of distance, we define the following functional
where is the identity on . Note that only represents the surface of one single vesicle .
3.3.2 The Model 3 and the Functional
Well-prepared with the above functional , we can now extend the single vesicle case - Model 2 to the multiple case - Model 3. Let for some . For each , we do the following problem:
where are the weight coefficients for and respectively. To translate this optimization problem: in fact, we are minimizing the Willmore energy , under the condition that, firstly, is hard to tough the region with a weight and secondly, is hard to get very closed to other vesicles with a weight .
The associated functional corresponding to this problem is given by
which is the augmented energy to be minimized in Model 3.
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 .
(Willmore Flow with Area Constraint and Obstacles - applied to Multiple Vesicles Case, Weak Form)
Suppose that are fixed as a weight coefficient for the functional and respectively. Given an initial shape , find the multiplier and the function according to the family of surfaces , such that on the time interval ,
for all test function .
The formulas of the shape derivatives of , and are clearly explained in the previous sections, so I only demonstrate the calculation of in this section. Since our model is implemented by FEM at the end, I estimate the shape derivative of in a tricky way. First, for any point , consider an open ball centered at with radius and define a function by
Since is independent on the geometry, we can apply Lemma 2.2.2 and then obtain the shape derivative of .
Applying this formula, then we can approximate the shape derivative of by piecewise implementation in FEM.
4.1 Time Discretization and Equation Split
4.1.1 Time Discretization
The model is implemented on the time-interval , though the final time can be chosen dependent on the stopping criteria set in the algorithm. Let
where , be a partition of the interval .
About Time Adaptivity: If one wants to obtain a more effective algorithm to reach the optimization shape of faster, it is more reasonable to make the time step
adaptive to the mesh size, because this FEM is using a moving-mesh. For me, I simply choose a comparatively small time step , which is fixed, for convenience. However, there is actually a goodness of using a small time step , 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 when is given. Recall the function defined by equation (3.4). From now on, we denote the numerical solution by at each time , which is viewed as the image of
with given approximation .
From equation (3.5) and (3.7), one can see that, if we solve directly, the order of the differential equation is high. To reduce this order, many previous work chose to split the formula . In our work, we solve a pair of unknown first, and then update . is defined as
and hence the following weak formula
for all test function . The function space of the test function will be discussed in the coming section.
4.2 Finite Elements
The following are some notations and definitions used in this section.
is the space containing .
is the parametrization space.
is a polyhedral approximation of , where is the triangulation of and the vertices of lie on .
is a -simplex in with its vertices .
is the vertex set that we use in FEM, which includes the vertices of and the mid-points of each edge of .
is the image of a function defined on such that is a polynomial of degree .
Denote . Then we have .
The function values of indicate the position of the vertices . The values of is obtained by following the rules:
If (i.e. is a vertex of the -simplex ), then . Since are on , then are also on .
If (i,e, is a mid-point of an edge of ), then is an orthogonal projection onto . Hence, also lies on .
Figure 4.2: Adjustment of the vertices
To make sure that the set is the orthogonal projection of on to , we adjust the position of the vertices whenever we obtain a new set from by solving the numerical problem.
By the above construction, one can conclude that all the points in are on . Besides, the existence and uniqueness of the polynomial satisfying the above two rules are proved by [16, Theorem 2.2.1].
, a -simplex in . We define the standard reference as the convex envelope with vertices .
In , .
In , .
represents the vertex set of .
For each -simplex in mentioned above, there exists a bijective mapping such that maps to .
Finite Element Space
The finite element space defined over the set is defined as
At each time step , is known. Hence, the finite element function space that we consider at the time step to find is .
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 for . After time discretization, we aim to solve (equation (4.1)) at each time . Instead of solving directly from , as discussed in Section 4.1.2, we split the formula by adding the equation (4.4), so we can solve first and then update by equation (4.2). The discretization (Section 4.2) gives us the space to solve . Based on these ideas, we can formulate the problems in Chapter 3 into the following discrete forms.
(Discrete Form of Problem 3.1.1.)
Suppose that are fixed as a weight coefficient for the functional and respectively. Given an initial shape , find and such that for each , ,
At each time step, is updated by
(Discrete Form of Problem 3.3.1.)
Given an initial shape , find and such that for each , ,
where the formulations of and 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 to be zero.
4.3.1 Linearization of and in FEM
Recall the formula for (3.14):
The function with smoothness (by applying exponential functions) is nonlinear. So is . In FEM, we use their linear approximation. They are simply made in the standard way:
Here is the Fréchet derivative operator. It is well-known that the approximation is good for if it is close enough to , so it is very natural that when we take the points from . When is small, the approximation could be good. Hence, we now replace by and by . More specifically, we write the approximation as: ,
Clearly, in the discretized form, we can replace by and then obtain the following semi-implicit formula for :