Remote-Url: https://math.nyu.edu/media/math/filer_public/4c/6d/4c6dcf27-3cce-44b7-9b3d-6365448f2286/tianruixu.pdf
Retrieved-at: 2023-05-27 22:55:18.390324+00:00

Mathematical Models for Red Blood CellTianrui Xu* and
Advisor: Charles S. Peskin*Applied Math Summer Undergraduate Research Experience (AM-SURE), Courant Institute of Mathematical
Sciences, New York University, New York, New York, United States of America
*tx392@nyu.eduAbstractThe red blood cell has a unique biconcave shape because this shape minimizes the bending
energy of the membrane. In this study, we want to find the simplest discrete model to generate
this shape. We use a triangulated surface to approximate a smoothly continuous surface and
re-define curvature on our discrete surface. We formulate an evolution equation involving
bending energy E, volume V , and area A, and use it to produce a biconcave shape that
simulates red blood cell.1IntroductionThe cell membrane is a very important part of the cell. It not only protects the cell from extreme
environments and virus invasions, but also connects the cell with the outer world. Nutrients are
transported into the cell through the membrane and waste is transported out. These transports
happen through different channels on a cell’s membrane and the exocytosis and endocytosis pro-
cesses. Armed with such functions, a cell’s membrane helps the cell to communicate with the outer
environment, although this communication is done without words [1].Due to the importance of the membrane, we would like to present a mathematical model that
can generate a cell’s membrane with the right shape. Among all types of cells, the red blood cell
has one of the most interesting shapes, which is biconcave. This special shape motivates us to build
a model for the red blood cell to show how the biconcave shape naturally emerges. Although there
are many past studies on models of red blood cell, the goal of this project is to find the simplest
model, and to use discrete differential geometry to generate a biconcave cell membrane from a
perfect sphere. This model can be further generalized to other membranes and surfaces that have
unique shapes.2 Mathematical Modeling2.1Initial ConditionTo get a biconcave cell, we start from a sphere. Then, we gradually pump out volume and
simultaneously let the bending energy to alter the initial shape such that the bending energy can
be minimized on the surface. Our model reversely mimics the experiment of putting the red blood1Figure 1(a)Figure 1(b)Figure 2(a)Figure 2(b)Figure 1(a)-(b) are examples of using triangulated surface to model a sphere. (a) has a total of 42 points and 80
faces and (b) has a total of 1280 points and 642 faces. The more flat triangles we use, the more smooth the surface
is. Figure 2(a) shows the three vertices with their coordinates Xf1 , Xf2 , and Xf3, along with the tetrahedron Tf
associated with the face f . Figure 2(b) shows the six tetrahedrons Tf1, Tf2, Tf3 , Tf4 , Tf5 , and Tf6 around the vertex k,
associated with the faces f1, f2, f3, f4, f5, and f6.cell in dilute solutions. When placed in dilute solutions, the biconcave red blood cell will become
spherical [6]. Therefore, transforming the cell from a sphere to a biconcave shape is experimentally
feasible. To get the intial shape of our model, we discretize the sphere as a faceted triangulated
surface. The more triangles we use, the more smooth the surface is [2], see figures 1(a) and 1(b).
After triangulating the surface, we track the coordinates of the points which are vertices of the
triangles.We use 642 vertices and 1280 faces. We give a unique index to every point and every face. We
store a point by assigning its coordinates to its unique index, and we store a face by assigning the
indexes of the three vertices of this triangle to the index of the face. In this work, Xk refers to the
Cartesian coordinates of the k-th point.2.2 Model SetupAfter initiating the cell membrane as a sphere, we want it to bend to a biconcave shape.
According to Canham’s discovery in 1970 [4], the reason a red blood cell is in a biconcave shape is
that the biconcave shape requires the least amount of bending energy. If we allow the membrane
of a cell to move freely, it will eventually bend to the biconcave shape to minimize its membrane’s
bending energy. Therefore, to obtain the shape of a red blood cell, we want our initial model,
the sphere, to change its shape over time, and during this process, we reduce the volume V and
minimize the bending energy E, while keeping its surface area A constant. To minimize the bending
energy E, we want to minimizeQ =dE
dt+1
2(cid:88)kdXk
dtdXk
dt,where t refers to time. Minimizing dE/dt gives us the direction that decreases E as quickly as
possible. To get a finite velocity we penalize speed by adding the second term (cid:80)
dXk
dt on the
right hand side of the equation for Q.dXk
dtkIn addition, we want to impose area and volume constraints. To preserve basic property of the
cell membrane, we want to keep the surface area A constant. This gives us dA/dt = 0. Given a
fixed surface area, the shape with the largest volume is a sphere. Since we start from a sphere
and want to keep surface area constant, we have to force volume V to decrease. This gives us the2second constraint dV /dt = dV0/dt, where V0(t) = V min +(V max −V min)e−at is a given function
with negative slope (V max refers to the initial volume, V min refers to the final volume, and a is
a positive constant).In summary, we would like to minimize Q, subject to dV /dt = dV0/dt and dA/dt = 0. Solvingthis problem, we translate it into the following system of ordinary differential equations (ODEs):dXk
dt= −− λ− µsubject to∂E
∂Xk∂V
∂Xk∂A
∂Xk,dV
dt=dV0
dtdA
dtand= 0,(1)where λ and µ are Lagrange multipliers for the constraints.2.3 CalculationAfter setting up the ODE equation, we would like to show how to calculate the volume V
enclosed by the faceted surface, the area A of the faceted surface, the bending energy E, and the
multipliers λ and µ. The 1280-sided polyhedron we construct to approximate a sphere is partitioned
into tetrahedrons. For each tetrahedron Tf (see figure 2(a)), its four vertices are the three vertices
of the triangle tf and the origin O. Face tf represents the f -th triangle face on the membrane’s
surface. Then the total volume V = (cid:80)
f Vf , where Vf = 1
6Xf1 · (Xf2 × Xf3) is the volume of the
tetrahedron Tf . Note that Vf is a signed quantity, and it is important that the vertices of each
triangle are numbered in counterclockwise order when viewed from outside the surface. This makes
Vf positive if the origin is within the region enclosed by the surface, but the formula is correct in
any case. When the origin is outside of the enclosed region, some of the Vf are positive and some
are negative and their sum is still the correct enclosed volume.The total surface area A is the sum of areas of all surface triangles and can be written as
f af , where af represents the area of the f -th triangle. We use Heron’s formula af =A = (cid:80)
(cid:112)s(s − a)(s − b)(s − c), where s = (a + b + c)/2, see figure 2(a).Bending energy E is defined as E = (cid:82) H 2da on a smooth surface, where H and a represent mean
curvature and area, respectively. By first variation of area, we can interpret mean curvature H of
a surface as the rate of change of its area, per unit area, when this surface is moving in the normal
direction to itself [5]. Therefore, we would like to apply this definition on our discrete surface to
define hk as the mean curvature of the vertex k, see figure 2(b). Then, E = (cid:82) H 2da ≈ (cid:80)
kak [3],
where k refers to index of vertices. We now focus on defining surface area ak and mean curvature
hk for the k-th vertex. Vertex k’s area is defined as the third of the sum of areas of faces around
vertex k. In other words, vertex k ’s area ak = (cid:80)
i∈F (k) (cid:107)A(i)(cid:107) /3, where F (k) is the set of indices
of faces that meet at the k-th point and A(i) is the surface normal of the i-th triangle. For the
f -th face in figure 2(a),k h2A(f ) =(Xf1 − Xf3) × (Xf2 − Xf3)=(Xf1 × Xf2 + Xf2 × Xf3 + Xf3 × Xf1).1
2
1
2Then hk should be the gradient of vertex k’s area, divided by its area, when vertex k is moving in the
direction that maximizes the rate of change of vertex k’s area. In other words, hk = (cid:107)(dak/dXk)(cid:107) ·
a−1
k . Then we are finally able to calculate E.The lagrange multipliers λ and µ are functions of time that ensure that the volume and areaconstraints in equation 1 are followed. To calculate λ and µ, we write the two constraints as3dV
dXk·dXk
dtdV
dt(cid:88)k
(cid:88)=== −k
dV0
dt,=dA
dt(cid:88)k
(cid:88)==dA
dXkdA
dXk··dXk
dt(cid:18)k(cid:88)kdA
dXk= −= 0.(cid:18)·dV
dXk
dV
dXk−∂E
∂Xk
∂E
∂Xk·k(cid:88)− λ− µ∂V
∂Xk(cid:88)− λkdV
dXk·(cid:19)∂A
∂Xk
∂V
∂Xk(cid:88)− µkdV
dXk·∂A
∂Xk−∂E
∂Xk
∂E
∂Xk·− λ− µ∂V
∂Xk(cid:88)− λkdA
dXk·(cid:19)∂A
∂Xk
∂V
∂Xk(cid:88)− µkdA
dXk·∂A
∂XkThis gives the system of equations(cid:80)
kdV
dXk· ∂V
∂Xk(cid:80)
kdV
dXk· ∂A
∂Xk− dV0dt − (cid:80)kdV
dXk· ∂E
∂Xk(cid:80)
kdA
dXk· ∂V
∂Xk(cid:80)
kdA
dXk· ∂A
∂Xk− (cid:80)
kdA
dXk· ∂E
∂Xk(cid:21)(cid:20)λ
µ ·= .2.4 RefinementSince we are using a discrete method to study a continuous problem, it becomes important
to keep the surface smooth. Smoothness is important to our model, not only because irregular
surface can produce errors, but also because the cell membrane has shear resistance. Therefore, to
restrict changes in the shape of the triangular facets of our model, we impose a penalty to our ODE
equation. Instead of solely minimizing E, we minimize E + α (cid:80)# of edges
l2
e, where le represents the
length of the i-th edge and α is the penalty’s weight. The weight α is chosen to balance (cid:80)
e l2
e and
√
√
E. Suppose that each triangle is equilateral, then (cid:80)
3A,
3 ≈ 2
√
where ne is the number of edges and nt is the number of faces. Then α should be 1/(2
3). By
minimizing E + (cid:80)
3, we limit the possibility that one edge will get extremely long compared
to other edges, and therefore limit the possibility that triangles will get deformed. Then our ODE
equation 1 becomes3 = 4af · 3nt/2e ≈ 4af · ne/e/2e l2e l2√√e=1dXk
dt= −∂E
∂Xk−1
√
3
2∂ ((cid:80)e l2
e)
∂Xk− λ− µ∂V
∂Xk∂A
∂Xk.(2)42.5 Solving ODEsWe use Euler method to solve equation 2:X (n+1)
k= X (n)k + ∆t ·dX (n)
k
dt,where ∆t is the time step, and X (n)
is the coordinate vector of the k-th point at time t = n · ∆t.
We apply equation 2 to calculate dX (n)
k /dt. The parameter ∆t is chosen small enough to avoid
numerical instability and to ensure that both volume and area constraints are followed. We use
∆t = 10−5 in our model.k3 Results5These nine figures are viewed from the side.Figures 3(a)-(i)Figures 3(a)-(i) show results at time t = 0, t = 1.5, t = 4.5, t = 7.5, t = 11.5, t = 15, t = 18, t =
20.5, and t = 22, respectively. We can see that the shape is changing gradually from a perfect
sphere to biconcave. It is important to note that the symmetry of the initial shape is broken in the
process of finding a direction to bend inward. Also, our refinement method is effective in keeping
the surface smooth at every time step. We can take a closer look at the final stage by looking at
figures 4(a)-(c). From figure 4(a), we can see that the surface is perfectly round viewing from the
top, which is exactly how red blood cell looks from the same point of view. Viewing from the side,
as shown in figures 4(b)-(c), the upper and bottom surfaces have bent inward as expected.Figure 4(a)Figure 4(b)Figure 4(a) is viewed from the top. Figures 4(b)-(c) are viewed from the side.Figure 4(c)Also, the model performs well in following constraints in equations 1 and 2. We plot the total
surface area A against time t, see figure 5(a). The plot shows that A has not changed much in
the time range [0 22], so the constraint on area dA/dt = 0 is followed. We also plot the total
volume V against time t, with the V0 curve defined in 1, see figure 5(b). Although two curves6are plotted, they overlap each other so we can only see one curve. This shows that V is changing
with V0 curve perfectly in the time range [0 22], so the constraint on volume dV /dt = dV0/dt is
followed. Figure 5(c) shows the distribution of side length ratios of triangles on the surface. Here,
we take each triangle and calculate the ratio of its maximum side length to its minimum side length.
The concentration of ratios around 1.0 shows that most triangles on the surface are approximately
equilateral; thus, deformations of triangles are restricted.Figure 5(a)Figure 5(b)Figure 5(c)Figure 5(a) shows the change of total surface area A over time t. Figure 5(b) shows the V0 curve (pink line) and the
change of total volume V over time t (blue line). Figure 5(c) shows the histogram of number of sides against side
length ratio.4 ConclusionUsing our discrete model, we successfully generate a biconcave shape that matches our expecta-
tion. Although there have been previous studies on simulating the shape of red blood cell, our model7have advantages in simplicity and can be easily applied to other shapes to obtain final phases that
minimize the surface energy. The refinement method used in this model also has practical values.
By equiping with the refinement tool, a continuous surface can be represented by discrete points
which will maintain the smoothness of the surface in any course of movement.5 AcknowledgementI would like to show my gratitude to my mentor Professor Charles S. Peskin for his great guidance
and encouragement. It is his wisdom that shows me the light through difficulties. I would also like
to present my special thanks to Dr. Pejman Sanaei and Mr. Jason Kaye for their patient assitance
and critiques of the research project. I would also like to thank Professors Miranda Holmes-Cerfon
and Aleksandar Donev for leading and supporting the undergraduate research program. Finally, I
would like to thank NYU Math Department for funding this research.References[1] Reece, J. B., Urry, L. A., Cain, M. L. 1., Wasserman, S. A., Minorsky, P. V., Jackson, R., andCampbell, N. A. (2014). Campbell biology. Tenth edition. Boston: Pearson.[2] Charles S. Peskin. Matlab codes that generate a sphere.[3] Charles S. Peskin. Notes on defining bending energy on a discrete surface.[4] Canham, P.B. (1970). The minimum energy of bending as a possible explanation of the bi-
concave shape of the human red blood cell. Journal of Theoretical Biology, 26(1), 61-81.
doi:10.1016/s0022-5193(70)80032-7.[5] Sullivan, J.M. (2006). Curvatures of Smooth and Discrete Surfaces. Discrete Differential Geom-etry Oberwolfach Seminars, 175-188. doi:10.1007/978-3-7643-8621-4 9.[6] Joseph, F. Hoffman (2016). Biconcave shape of human red-blood-cell ghosts relies on density
differences between the rim and dimple of the ghost’s plasma membrane. Preceedings of the
National Academy of Sciences, 113(51). doi:10.1073/pnas.1615452113.8