Physical Simulation and Visualization of Cells

(1)

Physical Simulation and Visualisation of Cells

Examensarbete utfört i Informationskodning vid Tekniska högskolan i Linköping

av

Johannes Nilsson

LiTH-ISY-EX--08/4071--SE

(2)

(3)

Physical Simulation and Visualisation of Cells

Examensarbete utfört i Informationskodning vid Tekniska högskolan i Linköping

av

Johannes Nilsson

LiTH-ISY-EX--08/4071--SE

Handledare: Niclas Wadstr¨omer Examinator: Robert Forchheimer

(4)

(5)

Presentationsdatum 2008-05-26

Publiceringsdatum (elektronisk version) 2008-06-17

Institution och avdelning Institutionen för systemteknik Department of Electrical Engineering

URL för elektronisk version http://www.ep.liu.se

Publikationens titel

Physical Simulation and Visualization of Cells

Författare Johannes Nilsson

Sammanfattning Abstract

Multicellular organisms develop from a single cell by repeated cell division and cell differentiation. This report presents a computer program for simulating dividing cells. The program simulates a single cell that goes through a series of divisions to create a small organism.

The evolving organism can be visualized both in a 3D view and in a lineage tree structure that shows the cell heradity. Information about both final cells and intermediary cells is easily available for the user.

A method for fast physical simulation of soft bodies based on a mass-spring system combined with a pressure model is used to simulate the physical properties of the cells. Fast collision detection for soft bodies is achieved using

hierarchical tree structure of axis aligned bounding boxes combined with an implicit penetration depth definition.

Biological control of cell division and differentiation is simulated using a model where the state of each cell is represented by

as a vector of substances present in the cell. The cell state is updated at discrete time steps by a model of DNA transcription.

Nyckelord

Physical, simulation, ,visualization, cell Språk

Svenska

X Annat (ange nedan)

Engelska Antal sidor 54 Typ av publikation Licentiatavhandling X Examensarbete C-uppsats D-uppsats Rapport

Annat (ange nedan)

ISBN -

ISRN LiTH-ISY-EX—08/4071--SE Serietitel -

(6)

(7)

Abstract

Multicellular organisms develop from a single cell by repeated cell di-vision and cell differentiation. This report presents a computer program for simulating dividing cells. The program simulates a single cell that goes through a series of divisions to create a small organism.

The evolving organism can be visualized both in a 3D view and in a lineage tree structure that shows the cell heradity. Information about both final cells and intermediary cells is easily available for the user.

A method for fast physical simulation of soft bodies based on a mass-spring system combined with a pressure model is used to simulate the phys-ical properties of the cells. Fast collision detection for soft bodies is achieved using hierarchical tree structure of axis aligned bounding boxes combined with an implicit penetration depth definition.

Biological control of cell division and differentiation is simulated using a model where the state of each cell is represented by as a vector of substances present in the cell. The cell state is updated at discrete time steps by a model of DNA transcription.

(8)

(9)

List of Terms

and Abbrevations

Term Description

AABB Axis aligned bounding box.

Autocrine signaling Self signaling, a cell is both signaling and receiving its own signal.

CellSim 3D The computer program described in this report that simulates dividing cells.

Collision detection The process of finding collisions between objects. Collision response The process of resolving collisions.

Contact dependent signaling Signaling between cells that are in direct mem-brane to memmem-brane contact.

DNA Deoxyribonucleic acid - stores genetic information in living organisms.

Endocrine signaling Signaling between cells by releasing hormones into the bloodstream that are later received by target cells far away.

Feature based methods Methods for collision detection that use features in objects to find colliding objects.

FEM Finite elements method. A method for solving dif-ferential equations over an irregular grid.

Gap junction A small channel between two cells that allow some particles to pass between them.

Gene A sequence of DNA that is necessary for generat-ing a functional product.

Hierarchical methods Methods for collision detection that use hierarchi-cal data structures to find collisions between ob-jects.

Implicit model A model where objects are defined by implicit functions.

Isopotential surface A surface defined as the points where a function has a certain value.

(12)

iv TERMS AND ABBREVATIONS

Term Description

Lineage tree A structure for showing the ancestry of cells. Lipid bilayer A double layer of lipid molecules that assemble tail

to tail. Such layers form the cell membranes. Mass-spring system A system where particles with mass connected by

springs are simulated.

OBB Oriented bounding box.

Paracrine signaling A form of signaling where the target cell or cells are close to the signaling cell.

Particle system A system where particles, possibly connected by constraints, are simulated.

Rigid body An object that doesn’t deform.

RNA Ribonucleic acid.

Shape matching A method for simulating soft bodies.

Soft body An object that can be deformed in response to external forces.

Spatial hashing A method for fast collision detection using a hash table.

Sphere tree A hierarchical structure used for collision detec-tion.

Synaptic signaling Fast signaling over long or short distances via synapses.

(13)

Chapter 1

Introduction

Multicellular organisms evolve from one single cell, through repeated cell division. As a result of differentiation at cell division and the environmental influence on the cells, cells will specialize themselves to perform different tasks in the organism. An adult organism contains many different kinds of cell types.

Methods that simulate the biological control of the cells and cell division are being developed. This master thesis will look at how this can be visualized in 3D and if this new representation of the cells gives new possibilities to add more parameters to the simulation.

1.1 Goal

The master thesis aims to use methods from computer graphics and computer games, and build a physical model of dividing cells. The output of the model can be used to visualize the simulation of dividing cells in 3D.

The main task is to find a suitable physical model to represent the cells. It shall handle deformable objects and include collision detection and collision response. Also the model shall be suitable for visualizing the results. Tests shall be made that shows how fast the choosen method is and how well it scales when the resolution of an object changes and when the number of objects in the simulation changes.

The result of the master thesis shall be a computer program where the physical model and the graphical visualisation of cells are combined with a model of the biological control of the cells. A short video or image sequence shall be created that shows the first cell divisions of a simple modelled organism.

The model can also give possibilities to introduce some new parameters and see how these affect the evolving modelled organism. It is good if the computer program is designed so that tests with different spatial parameters can be made

(14)

2 CHAPTER 1. INTRODUCTION

in an easy way.

Different strategies for placing the division plane when a cell divides is one of these parameters. The policy can be to divide a cell according to some globally defined axis or to let a cell divide in the direction where it has the largest spatial extent. One simple policy is to choose a random division plane for each division and that policy can serve as a reference.

Another parameter of interest is how the cells connect to each other. Two distinct ways of doing this can be tested, either to connect cells at division or to let a special gene expression signal to the cell to make connections to its neighbours.

1.2 Outline

Here a short outline of the six chapters in this report is given. Chapter 1 gives an introduction and presents the goal of the master thesis. Chapter 2 introduces cells and gives some basic knowledge about them needed to read the following chap-ters. Chapter 3 gives an introduction to soft body simulation and discusses pos-sible methods for physical simulation and collision detection. Chapter 4 presents the implemented method for simulating cells. Chapter 5 presents the resulting program for physical simulation of dividing cells and also shows the results from different test runs. In chapter 6 the results are discussed and some possible future work is proposed.

(15)

Chapter 2

Cells

This chapter will give a short introduction to what cells are and introduce some properties of cells that are of importance for this thesis. The study of cells is a huge scientific area that cannot be summarized here. To get a more complete introduction to cells see the books by Alberts[1] and Prescott[11].

2.1 Function and structure of cells

Cells are the smallest building blocks of all living organisms, both animals and plants. Some simple organisms like amoebas or bacterium consist of a single cell while a complex organism like an adult human being has about 1014

cells. A fundamental property of cells is that they have the ability to make copies of themselves by cell division [11].

Figure 2.1: A schematic image of a cell. The plasma membrane is the outer boundary of the cell. Inside the cell is the cytosol and the organelles. Most cells also have a nucleus.

A plasma membrane defines the outer limit of a cell and consists of a lipid bilayer. It separates the cytoplasm inside the cell from the extracellular space

(16)

4 CHAPTER 2. CELLS

between cells. The cytoplasm consists of a water-like solution called cytosol and the cytoplasmic organelles that can be seen as the organs of a cell. Cells can also have a nucleus where most of the DNA and RNA contents are created and stored. The nucleus has its own nucleic membrane that separates it from the cytoplasm [1].

2.2 Cell types

When studying cells they are often divided into prokaryotic cells and eukary-otic cells. Prokaryeukary-otic cells are less complex than the eukaryeukary-otic and contains less DNA, they don not have a nucleus. Bacterium are prokaryotic cells. All prokaryotic organisms are unicellular [11].

Multicellular organisms are built up by eukaryotic cells. In a multicellular organism there are many different cell types, like muscle cells, epithelial cells and nerve cells. Over 200 identified different cell types can exist in the same organism and they perform a wide variety of tasks [1].

2.3 The cell cycle and cell division

A cell reproduces by going through a series of events in which it duplicates its contents and then divides in two. This series of events from when the cell is created by duplication to when it divides again is called the cell cycle [1].

Figure 2.2: The four phases of the cell cycle. Cell division occurs in the M phase, G1, G2 and S are called the interphase. During the interphase a cell

roughly double its size and copies its DNA contents.

Figure2.3shows the phases of the cell cycle. In the M -phase the cell divides by nuclear division followed by cytoplasmic division. The G1, S and G2 phases

are together called the interphase. Cells in an adult organism roughly double their size during the interphase, but in the early development of an organism

(17)

2.4. GENES AND GENE EXPRESSION ₅

cells will divide without growing. In the S-phase the cell duplicates its DNA contents [1].

2.3.1 Lineage tree

One way of looking at the development of an organism is to view the lineage tree. This shows how different cells are related to each other. In the tree only the leaf nodes are currently existing cells, while the internal nodes represent ancestor cells.

Figure 2.3: The lineage tree of the first cell divisions for the nematode C. elegans.

2.4 Genes and gene expression

Cells store their genetic information in DNA. DNA is a very long double-stranded molecule that stores hereditary information as a linear sequence of the four monomers A, T, C and G [1].

A gene is a complete sequence region in the DNA that is necessary for gen-erating a functional product. A functional product can be a protein or a piece of RNA [13].

The DNA code is practically identical for all cells in an organism while it is the expressed genes that differs between the different cell types [1]. Of all the genes that are encoded in the DNA of a cell, only some are expressed. Genes are expressed through the protein products or RNA molecules they encode. Expres-sion of a gene can be controlled by transcriptional activators and transcriptional repressors that bind to the binding sites for the gene [1].

(18)

6 CHAPTER 2. CELLS

2.5 Cell communication

Cells communicate in many different ways. Extracellular signaling is performed by secreting signal molecules into the extracellular space. These molecules then bind to specific receptors on the receiving cell, which in turn sets of a chain of events that results in some kind of response [1].

Figure 2.4: (A) Contact dependent signaling requires the cells to be in direct membrane to membrane contact with each other. (B) Paracrine signaling lets cells communicate locally by secreting substances into the extracellular space. (C) Synaptic signaling is performed by sending electrical signals along a process called an axon. At the end of the axon the electrical signal is converted back to a chemical signal. (D) Endocrine signals are passed through the blood stream of a larger animal.

Some signaling between cells is contact dependent because the signaling molecules are bound to the plasma membrane of the signaling cell. In that way only cells that are in direct membrane to membrane contact can be affected. Paracrine signaling is done by secreting signal substances into the extracel-lular space and let them diffuse there. The paracrine signaling is only local. Communication over longer distances is done by synaptic signaling where neu-rons have long processes called axons that can transmit electrical impulses long distances. At the end of an axon there is a neurotransmitter that converts the electrical impulse back into a chemical signal at a cell junction called a chemical synapse. Another way of signaling over long distances is endocrine signaling where hormones are released into the bloodstream by the signaling cell [1].

Many of the substances used in paracrine, endocrine and synaptic signaling are the same. The difference is the speed substances can be delivered to the target cell and how selective the signaling can be [1].

(19)

2.6. CAENORHABDITIS ELEGANS ₇

Autocrine signaling is when a cell sends a signal and also has receptors so it can receive the signal. This is often used in the development of an organism when groups of cells take group decisions about how to differentiate. If most of the cells in a neighborhood sends out the same autocrine signal that signal will be the strongest there and all the cells then know how to differentiate [1].

Cells that are in contact can create gap junctions that are small channels between the plasma membranes of the cells that allow some small particles to be shared by both cells. This might have similar effects as the autocrine signaling has in development decisions because it makes cells aware of the state of their neighbors [1].

2.6 Caenorhabditis elegans

Ceanorhabditis elegans is a small nematode worm and is about 1 mm long. The adult worm has about 1000 cells and its development is very well studied in terms of genetic control [1].

There are several reasons why C. elegans is beneficial to use when studying cell development. The worm is very well studied and there are lots of information available. In 1998 it became the first multicellular organism for which a complete genome sequence was obtained. The complete genome sequence is about 100 Mb and it is well annotated. It is also easy to generate mutations of C. elegans [5].

The fate of all cells in C. elegans is almost perfectly predictable. C. elegans life begins with one cell, the egg. Through repeated cell divisions the egg cell gives rise to a small worm consisting of 558 cells inside the egg shell [1].

(20)

(21)

Chapter 3

Soft body simulation

Today most computer game engines only implement rigid body physics and there exist several commercial software packages that efficiently simulates rigid bodies. Soft body simulation on the other hand is more often used in non-real time simulations where methods like the finite elements method, FEM, can handle both plastic and elastic soft bodies. Recent research shows that soft body simulation in real time is possible on modern computers. These simulations can be used in areas like virtual surgery and computer games. New methods can in some cases be used to speed up simulations or make more complex simulations possible.

This chapter will give an introduction to simulation of soft bodies. There are several possible approaches and some will be discussed here. Johansson [7] compares different methods for simulating soft bodies in computer games and suggests a method based on shape matching. Matyka [10,9] proposes a model based on a mass-spring system combined with a pressure model to simulate soft bodies.

3.1 Modeling soft bodies

When considering which model of soft bodies to use it is a good thing to first look at what the application is. The finite elements method is quite general and can be used in a wide range of simulations while other more specialized models are suited for a specific simulation. A more specialized model often means increased performance of the simulation.

The following three sections will describe three groups of methods to simulate soft bodies. In particle systems a complex object is represented as a group of particles that are connected in some way to build up the object. Implicit models

(22)

10 CHAPTER 3. SOFT BODY SIMULATION

use continuous functions to represent objects. The finite elements method is a fast way of solving partial differential equations over an irregular grid in an efficiant way and can be used to simulate soft bodies.

3.2 Particle systems

In rigid body dynamics a number of parameters such as rotation, translation, velocity of the center of mass and angular velocity can describe the state of a large and sometimes complex object. This approach is often difficult when handling soft bodies because of the problem of updating and keeping track of changes in the object at run time. If the deformations are unconstrained it can be hard to define the state of the object. In this case it may be a better approach to break down the object into simpler parts and add some sort of connections that give the object the right properties.

A point particle with mass can be simulated very fast and the complexity of the simulation depends mostly on how the particles are connected. Equation3.1

is used to calculate the acceleration of each particle. It simply divides the force induced by the simulator by the mass of the particle.

ai =

fi

mi

(3.1)

3.2.1 Mass-spring systems

In mass-spring systems the particles that represent an object have a mass and are connected by springs. A damped spring can be described by equation 3.2, and the acceleration is calculated in equation 3.3. In the equations ks is the spring

constant and kdis the damping constant.

am + kdv + ksx = 0 (3.2)

a = −ksx

m −

kdv

m (3.3)

Mass-spring systems are often used in simulation of cloth and other 2D ob-jects. When simulating 3D objects the problem of finding springs and coefficients that gives the object the desired physical properties is hard. Also the number of springs needed to simulate an object increases rapidly when the dimension in-creases. Going from two to three dimensions in the case of a square and a cube, only six springs are needed to connect a square symmetrically compared to the 28 springs that are needed to connect a cube symmetrically.

A mass-spring system is only conditionally stable. That means that depend-ing on the time step used in the simulation and the method for updatdepend-ing the differential equations3.2and 3.3the system may be stable or unstable. Different

(23)

3.2. PARTICLE SYSTEMS ₁₁

methods can be used to avoid instable systems, more about that is discussed in section 3.5.

Figure 3.1: Particles connected by springs.

3.2.2 Shape matching

Fundamental to shape matching is that objects have a base shape that is defined by the position of a set of particles when the object is at rest. Each particle in the simulated object has a corresponding particle in the base shape.

Figure 3.2: In shape matching every point in the object (green dots) gets a dis-placement vector that is directed toward the corresponding point in the translated and rotated resting shape (red dots).

When simulating a soft object with the shape matching method physics are first applied to the particles as if they were independent. Then the best matching translation and rotation of the base shape compared to the current positions of the particles is calculated. For each particle a penalty force is added to move the particle towards the corresponding particle in the rotated and translated base

(24)

shape. The force may be a linear force proportional to the distance, a damped spring force or any other force that will make the particles converge towards their position in the translated and rotated base shape.

To calculate the best matching translation and rotation of the base shape the first step is to translate the base shape so that the centers off masses are the same. The best rotation is defined by the rotation matrix that minimizes the sum in 3.4, where xi are the points in the deformed object, yi are the corresponding

points in the base shape and mi are point masses.

X

i

mi(R(xi− xcm) − (yi− ycm))2 (3.4)

Finding the rotation matrix R is difficult, but if the problem is relaxed to3.5

where A is a linear transformation matrix the problem is solvable. An approxi-mative rotation matrix can then be found by RS-decomposition of the transfor-mation matrix. How this is done is described by Johansson [7].

X

i

mi(A(xi− xcm) − (yi− ycm))2 (3.5)

Control of volume and placticity of the object is also possible by utilizing the information in the transformation matrix A. This is also described by Johansson [7].

3.2.3 A pressure model for soft bodies

The pressure model for soft bodies is proposed by Matyka [10]. It addresses the problem when using mass-spring models for 3D objects. The idea is that a 2D shell of springs is combined with a model of inner pressure in an object. This way no inner springs has to be used. The pressure can be calculated using the ideal gas law 3.6.

P V = nRT (3.6)

In the equation P is the pressure, V the volume, n is the mol number of the gas, R is the ideal gas constant and T the temperature. Under the assumption that nRT is constant the pressure can be expressed in equation 3.7.

P = kp

V (3.7)

The force induced on a face by the pressure, given by equation4.3is propor-tional to the area of the face, Ai, and directed in the surfaces normal direction,

ˆ ni.

~

(25)

3.3. IMPLICIT MODELING ₁₃

Figure 3.3: The inner pressure adds extra forces on each vertex.

Every point in an object receives pressure forces from all its neighbouring faces. The pressure model approximates a thin elastic shell filled with gas. By setting the pressure constant and the spring constants will control the properties of the object in a more intuitive way compared to pure mass-spring systems.

3.3 Implicit modeling

Implicit modeling does not depend on a mesh to define objects, instead continuous functions define objects. Often a sum or other combination of simple functions is used to describe an object. The isopotential surface, given by equation 3.9, of an object is defined as all points where a field function has a given value [3].

{x ∈ R3

|f (x) = c}, (3.9)

Equation3.10 defines the field function as a sum of functions fi(x).

f (x) =

n

X

i=1

fi(x). (3.10)

Outside the object f (x) > c and inside the object f (x) < c. This means that the value of the potential function in equation3.11can be used to determine how far from or how deep inside an object a point is. This is very useful for collision detection. By calculating intersections between implicit models collision forces can be determined.

p(x) = f (x) − c (3.11)

A result of using analytical field functions to define objets is that the surface has infinite resolution. To visualize the surface a raytracer can be used or the surface can be approximated by a mesh using methods like the marching cubes method [8].

(26)

3.4 Finite elements method

The finite elements method, FEM, is a general technique for creating approximate solutions to boundary-value problems. The method divides the domain of the solution into a finite number of elements and solves the equations approximately over the collection of elements [2].

FEM is based on the continuum mechanics, but it has been used in a wide range of scientific areas including physical simulations of deformable objects. It can be used to simulate both plastic and elastic objects and it can handle objects where the material properties varies over the object [7].

When simulating 3D objects using FEM tetrahedrons are most often used as finite elements, so the object that shall be simulated needs to be divided into a set of tetrahedrons. This can be hard to do without using too many tetrahedrons. The tetrahedral meshes on the other hand gives an advantage compared to triangle meshes when it comes to collision detection [7].

3.5 Differential equations

Differential equations in the physical simulator are integrated numerically. There are several possible ways of doing this. Euler integration first updates the posi-tions based on velocity and then updates the velocity based on the acceleration in the way showed in equation 3.12.

(

xt+1 = xt+ vt∆t

vt+1 = vt+ at∆t

(3.12)

Verlet integration use the two previous positions and acceleration to calculate the new position using equation 3.13. The Verlet integration method has better stability properties compared to Euler integration. Also another good feature with Verlet integration is that the current position can be moved by for instance a collision handling and velocities will change correspondingly, without increasing the energy in the system [6].

xt+1= 2xt− xt−1+ at(∆t)2 (3.13)

3.6 Collision handling

To handle collisions in a physical simulation, the collisions first need to be de-tected in a procedure called collision detection and then resolved in a procedure called collision response. Depending on how collision response is implemented the collision detection procedure may need to provide more information than just if objects are colliding or not. Here it is assumed that the collision

(27)

detec-3.6. COLLISION HANDLING ₁₅

tion mechanism provides information about penetration depth if objects are in a colliding state.

Handling collisions in a fast and correct manner is perhaps one of the most difficult tasks when dealing with deformable objects. One reason that methods developed for rigid bodies cannot be applied to deformable objects is that those methods are based on the assumption that the object is convex or can be divided into convex parts before the simulation starts. This is generally not possible if the objects are deformable.

Methods for collision detection can be divided into subgroups. Feature based methods take advantage of coherence in time in a simulation. That means that there is a little change in position of objects between two frames and thus a feature that separates two objects in one frame will be likely to do so in the next frame too. Hierarchical methods use trees of bounding objects to break down the problem of collision detection so that objects that are not in a colliding state can be efficiently removed from a list of possibly colliding objects.

3.6.1 Feature based methods

One feature based collision detection method is the separating plane method. For a pair of convex polygon meshes that are not in a colliding state there always exists a separating plane for which all points in the first mesh are on one side and all points in the other object is on the other side. One way of finding such a plane is to test all faces on one of the objects to see if the plane coresponding to the face separating plane. Once such a plane is found it is stored in the object and the next frame begins with testing if that plane separates the objects again. If there is a great time coherence it most often will be a separating plane again. That is what makes the method fast.

This only works for convex polygonal models or models that can be divided into convex parts. For rigid objects the slow process of dividing the object into convex subparts can be done before the simulation, but for deformable objects it must be done at run time.

3.6.2 Hierarchical methods

Hierarchical methods for collision detection store the physical data of an object in a hierarchical structure that is fast to search for collisions. The basic idea is that the first level of collision detection tests a simple bounding volumes that contains the whole object, if that test shows no collision the search can be ended directly. If the test shows that the bounding volumes collide the search continues until either there is no collision or the finest level of the collision search is reached. Common bounding volumes are bounding spheres, axis aligned bounding boxes , AABBs, and oriented bounding boxes, OBBs.

(28)

Sorkin [12] describes hirerarchical methods for collision detection and presents a method using sphere trees combined with a good heuristic for updating the tree after deformation which is important for objects that has a lot of deformation in at each time step. The principle is to begin the update at the leaf nodes and traverse the tree upwards to avoid multiple calculations at each node.

Gino van den Bergen [15] proposes a method based on AABB trees for efficient collision detection of deformable objects and states that if lots of deformation occur each frame AABB trees is a better choice of method compared to OBB trees. This is because updating an OBB is a more complex operation than updating an AABB.

3.6.3 Spatial hashing

Teschner and Heidelberger [14, 4] describes methods for consistent penetration depth estimation and collision response for flexible objects called spatial hash-ing. The method divides space into discrete regions with a specified size and an index in each dimension ix, iy and iz. These indices are used to calculate a hash

value that corresponds to a 1D hash table, see figure 3.4.

Figure 3.4: The spatial hashing method maps volumes in 3D-space onto a 1D hash table.

The first pass of the collision detection algorithm inserts each vertex into the hash table and the second pass tests each primitive against the hash table to find possible collisions. If a hash collision is detected a real collision test is performed. An advantage of this method is that self-collisions can be computed in a fast way.

3.6.4 Implicit methods

For implicit models the field function can be used directly to determine if a point is inside an object and the gradient of the field function can be used as a direction for collision response. This makes collision of points against an implicit model very efficient.

(29)

3.6. COLLISION HANDLING ₁₇

3.6.5 Collision response

When a collision is detected and the penetration depth calculated, some action needs to be taken to move the colliding objects apart. A common way to do this is to add a penalty force to all the penetrating particles that moves them out of the other object the next time they are updated. If Verlet integration is used it is possible to directly move the particles out of the other object without loosing any energy in the system. In both cases friction may be added as a force that slows the velocity of the particle in the plane perpendicular to the penetration vector.

(30)

(31)

Chapter 4

Method

This chapter describes the implemented method for simulating cells. First an overview of the simulation process is given and then each step is explained in more detail.

4.1 Overview of the simulation

The simulator is divided into a physical simulation unit and a simulation unit for control of cell division and differentiation. Each unit with its own defined interfaces as shown in figure 4.1. These are controlled by a simulator unit that manages the cell objects. The simulator unit also keeps track of the cell lineage tree as the cells divide.

Figure 4.1: Structure of the simulator.

The process of updating cells is controlled by the simulator unit. Collision

(32)

20 CHAPTER 4. METHOD

handling is handled by a separate unit that will update objects with penalty forces to resolve collisions. Physical simulation requires a much smaller time step compared to the biological simulation. This is handled by the method wantsT oRunCellcycle() in the biological cell unit. Figure4.2outlines the process of updating cells.

Figure 4.2: The steps that are performed to update a cell object.

4.1.1 Interfaces

The cell object is built out of a biological cell unit and a physical cell unit as shown in figure4.1. An interface has been defined for communication between the physical and the biological cell. The biological cell unit shall implement the meth-ods described in table 4.1.1. The methods wantsT oDivide() and wantsT oDie() returns a boolean value if the cell wants to divide or die and leaves to the simula-tor to handle that action. The method wantsT oGrow() returns a positive value if the cell wants to grow and a negative value if the cell wants to shrink, otherwise a zero value is returned. The method wantsT oConnect() returns a positive value to indicate that it wants to connect and a negative value if it does not want to connect. It is the combined value of two cells that shall be used by the simulator to decide if a connection shall be made.

Method name Description

wantsToDivide() Returns true if the cell wants to divide.

wantsToDie() Returns true if the cell wants to die (programmed cell death). wantsToGrow() Returns a value if the cell wants to grow or shrink.

wantsToConnect() Returns a value that indicates if the cell wants to connect. Table 4.1: Interface for biological cell unit.

The physical cell unit implements the methods described in table4.1.1. This interface defines the parameters available for the biological cell unit to make

(33)

4.2. SOFT BODY SIMULATION ₂₁

decisions based on physical properties.

Method name Description

getMass() Returns the cell mass. getVolume() Returns the cell volume.

getPressure() Returns the pressure inside the cell. Table 4.2: Interface for physical cell unit.

4.2 Soft body simulation

The implemented soft body simulator is based on a mass-spring system with a pressure model added. That is both an efficient model computationally and has some physical relevance to the cell as it simulates a soft shell filled with a fluid. A triangulated sphere, see figure4.3is used as the mesh that defines the deformable cell object. The sphere object can be replaced by any other triangle model, but it is assumed that the model is a closed, acyclic mesh with approximately equal sized triangles.

Figure 4.3: Wire frame view of a cell.

Figure 4.4shows the steps that are performed for each time step in the soft body simulator. An extra spatial constraint is added so that the cells are forced into the unit sphere. This simulates a spherical egg shell in which the organism grow.

4.2.1 Mass-spring system

The mass for each particle is calculated as the total mass of the object divided equally to each particle. This is only a good approximation if the vertices are

(34)

Figure 4.4: Outline of a time step in the physical cell simulator.

evenly spread in space.

mp =

mtot

Nvertices

(4.1)

Each edge in the mesh is considered as a spring with the rest length set to the length of the edge when the model is loaded. The spring constants ks and kd

controls the stiffness of the cell membrane. Equation 4.2gives the force induced by a spring, ∆l is the difference in length of the spring and ∆v is the relative velocity of the particles in the direction of the spring. The force is directed in the direction of the spring.

fs= ks∆l + kd∆v (4.2)

4.2.2 Pressure force

The pressure model proposed by Matyka [9] is used with the assumption that there is also a pressure outside the cells so that the difference in pressure is near zero when a cell is at rest. This is given by equation 4.3 where P0 is the outer

pressure

P = kp

V − P0. (4.3)

The volume of a mesh object can be calculated using equation4.4, appendix

A shows details about this. In the equation vi are vertex positions and Fi a face

(35)

4.3. COLLISION HANDLING ₂₃

V =X

iǫS

(v1+ v2+ v3)Fi

6 · n(Fi) × A(Fi) (4.4)

The pressure force is distributed to each face and proportional to the face area. For each particle equation 4.5gives the force induced to that particle from the adjacent faces, Fi is a face adjacent to the point i and AFi is the area of the

face Fi. fpi= X Fi P AFi 3 (4.5)

4.2.3 Updating the equations

The verlet integration scheme has been used because of its stability properties and because it allows that constraints are resolved by simply moving particles. To be able to resolve the spring equations the velocity also needs to be calculated at each time step. This is performed by simply taking the difference between the last and current position. Equation 4.6 gives the new position and equation4.7

calculates the velocity for a particle.

xt+1= 2xt− xt−1+ at(∆t)2 (4.6)

vt+1= xt+1− xt (4.7)

4.2.4 Scaling objects

Scaling of objects can both mean changing the volume of the object and changing the length of the springs and thus change the area of the surface of an object.

In the case of cell division the volume is divided into two equal parts. If nothing is produced at cell division both the sum of volumes and the sum of the surface area of the two daughter cells should be equal to the parent cells surface area and volume. In practice that would mean that the plasma membrane would be stretched tighter and tighter at each division. It is realistic to assume that new membrane is produced when cells divide to keep the pressure inside a cell constant and near zero. Alternatively a cell could loose volume at division to avoid increasing the inner pressure. The implemented method adds extra cell membrane with a given factor to prevent increasing pressure.

4.3 Collision handling

Collision handling is performed in two steps. First collision detection is done which means that the penetration depth of each vertex is calculated. Second, collision response applies forces to resolve the colliding state. The separation of

(36)

these tasks is important to avoid getting a simulation that is dependent on the order in which the objects are traversed.

4.3.1 Collision detection

The implemented method for collision detection is based on a definition of pene-tration depth that use implicit functions combined with an AABB-tree for each object. Each vertex has a local influence that is approximately the length of an edge, which means that a point inside a triangle is affected by all three corner points. The penetration depth given by each point is the distance to the vertex normal plane with sign. Figure 4.5 shows the resulting penetration depth for circles with different resolution. How much each point influence the penetration depth calculation depends on a weight function that is inversely proportional to the distance to the point.

Figure 4.5: Penetration depth for different resolutions of a circle. (A) With N = 3 gives a rough approximation of the surface. (B) Increasing the resolution of the model to N = 5 refines the approximation. (C) At a higher resolution N = 8 the surface of the object is smooth, but on the other hand the estimation of depth is inaccurate at the center of the object. This can lead to problems if objects collide at high speeds compared to the time step of the simulation, so a good choice of resolution must be made when using this method.

Because each vertex only has local influence an AABB-tree can be used as a hierarchical structure to search for collisions. This gives fast collision detec-tion and easy calculadetec-tion of penetradetec-tion depth. Figure 4.6shows the process of calculating the penetration depths.

4.3.2 Collision response

When a penetration depth is given by the collision detection algorithm, collision response is performed to resolve the collision state. The approach chosen is to add penalty forces that will move the vertices out of the other object in the next integration step. The penalty force is calculated by equation 4.8. This force can also be seen as a normal force and used to implement friction.

(37)

4.4. CELL SIMULATION ₂₅

Figure 4.6: Collision handling. (A) Two intersecting objects with vertices marked. (B) A closeup of two vertices that are intersecting and their bound-ing spheres and normals. Penetration vectors, marked with ’p’, is penetration depths calculated from each vertex in the other object. (C) The resulting pene-tration vector is the weighted average of the two penepene-tration vectors, the velocity and resulting penalty force is also shown.

fpenalty = −0.5 ˆd |d| (∆t)2 + v · ˆd ∆t + fother· ˆd ! (4.8)

4.4 Cell simulation

A unit for simulation of biological cell control has been implemented with the functionallity needed to produce dividing cells. This unit has been designed to make it possible to implement other different cell models later.

4.4.1 Modeling genes and DNA

The notation of logical genes proposed by Robert Forchheimer is used here. A gene gi(ak|al, am) is said to be gene gi and codes for transcription of substance

ak given that substances al and am is present, if that is the case the gene is said

to be active.

A modeled organism has a DNA code that is defined as a list of genes. The list of genes in the DNA is constant and equal for all cells in the same organism. The state of a single cell is described by the substances present. At each cell cycle

(38)

the list of genes is traversed. If all required promoting substances for a gene are present in the state vector of the cell it will transcribe its product that then will be inserted into the next state vector. A substance must be produced each time step to stay in the state vector.

Special substances initiates cell division, cell death and asymmetry in the daughter cells after cell division. Division and death substances must be produced by genes in the DNA code, while asymmetry substances are produced by the simulator into the state vector at each cell division.

4.4.2 Cell division

When a cell divides in the simulation its DNA content and state vector will be copied to both daughter cells. Also the simulator inserts special substances into each daughter cell to make them unique, in that way asymmetry is achieved.

The cell object implements a tree node interface and at divison the two new daughter cells are placed as children nodes to the parent cell node. Then the parent node is removed from the physical simulation but remains in the lineage tree so that it can be visualized.

Figure 4.7: Cell division: (A) A division is placed through the cell objects center of mass. (B) The cell object is copied and points in the new cell objects are copied and projected to different sides of the division plane.

Figure 4.7 shows how cell division is performed in the simulator. A division plane is created through the cell objects center of mass, then the cell object is copied and the points in the new cell objects are projected to different sides of the division plane.

(39)

4.4. CELL SIMULATION ₂₇

4.4.3 Cell connections

Cell connections in the simulator is implemented as adhering forces. For cells that are connected, the force is constant and directed towards the other cell. Equation4.9gives the force if ˆv is the unit vector directed towards the other cell.

fc = kcˆv (4.9)

The simulator holds a list of all the cell to cell connections. An important consideration when making connections is the policy for handling existing cell connections at cell division. One reasonable policy is to let the divion plane decide which of the existing connections that goes to one or the other of the daughter cells. Another possibility is to keep all connections for both daughter cells and instead let cells release connections under some conditions.

(40)

(41)

Chapter 5

Results

In this chapter the results from the physical cell simulation are presented. First the implemented computer program for simulating cells, CellSim 3D, and its functions are described and then results from running simulations are presented.

5.1 CellSim 3D

CellSim 3D is the resulting computer program for simulating cells. It can visu-alize both the lineage tree of the simulated cells and the physical appearance of the resulting organism. Figure 5.1 shows the main window of CellSim 3D. The program is implemented using Java with swing components for creating the user interface and OpenGL for the 3D rendering.

Figure 5.1: The main window of CellSim 3D. Left is a 3D view of the cells in the simulation and to the and to the right is the cell lineage tree.

(42)

30 CHAPTER 5. RESULTS

5.1.1 Physical simulator

Test runs with the physical simulator shows that it scales well up to about 128 cells, but above that the performance is poor. Timings are listed in table 5.1.1.

Cells Time (ms) Time (ms) / cell

8 14 1.75 16 30 1.88 32 62 1.94 64 116 1.81 128 241 1.88 256 6061 23.68

Table 5.1: Update times for different number of cells in the simulation. This shows that the simulator scales well for up to 128 cells.

5.1.2 Biological simulator

In the implemented model of cell control no interaction from the physical model will affect the resulting cell lineage. This makes it easy to generate code for a known lineage tree. The DNA code is specified in a text file where a line can specify a start vector, a gene or be a comment.

The following line specifies a start vector with a0, a1 and a2 present in the

first cell:

START: a0, a1, a2

This could also be written as three separate lines:

START: a0 START: a1 START: a2

To specify a gene g1(a3|a1, a2) the following syntax is used:

g1(a3|a1, a2)

Comments are lines that begin with a #:

# This line is a comment

5.1.3 Visualization

CellSim 3D has possibilities of visualizing both the physical cell simulation and the lineage tree structure produced by the cell divisions. There are three parts of the interface; the 3D view, the lineage tree view and the information box.

(43)

5.2. SIMULATIONS ₃₁

In the 3D view the mouse is used to rotate and zoom the virtual camera around the cells. Clicking on a cell in the 3D view will select it in the lineage tree and show information about it in the information box. It is possible to toggle the display of vertex normals and forces from the physical simulator in the 3D view and also to toggle between wire frame and solid view. This is useful to see cells inside a large organism.

The lineage tree view shows the cell heradity information. Clicking on a leaf node will select the cell in the 3D view and show information about it in the information box. Selecting an internal node in the tree will select all the cells in the 3D view that originate from this cell.

Figure 5.2: The lineage tree of C. elegans with cells originating from the cell EMS colored red, cells originating from P’ colored green and cells originating from AB colored blue. The color coding makes it easy to track where cells are located.

CellSim 3D has a method for cell coloring depending on substances in the state vector. A list of color associated substances is given as input and if a cell has more than one colored substance the colors are mixed. One way of obtaining this is to add an extra gene somewhere in the tree that produces a coloring substance to that subtree. Figure 5.2 shows an example of a tree where three different subtrees have been given different colors.

5.2 Simulations

This chapter presents the result from simulations with CellSim 3D. Two testruns with small cell trees and one with the lineage tree of C. elegans are performed.

(44)

32 CHAPTER 5. RESULTS

5.2.1 Simple tree with unique nodes

In this testrun a simple tree that generates four cells, shown in figure 5.3, is simulated. The resulting lineage tree and physical appearance are also shown in the figure.

Figure 5.3: A simulation where all the final cells are unique. The input tree, lineage tree from CellSim 3D and the 3D view of the final cells are displayed. Colors have been assigned to the substances to make the 3D view match the input tree.

The tree has four unique final nodes and three internal nodes, a total of 13 genes are used to code this tree. The complete DNA code is presented below:

# Start vector START: a0

# Internal node at level 1 g1(a1|a0)

g2(divide|a0)

# Internal nodes at level 2 g3(a2|a1,A)

g4(a3|a1,B) g5(divide|a1)

(45)

5.2. SIMULATIONS ₃₃

g12(a6|a6) g13(a7|a7)

5.2.2 Simple tree with two equal subtrees

In this test run a small tree, shown in figure 5.5is simulated. The tree has five final nodes and two equal subtrees. In the figure the input tree, resulting lineage tree and physical appearance are displayed.

Figure 5.4: Simulating a simple tree with two equal subtrees.

Figure 5.5: A simulation where two subtrees are equal. The input tree, lineage tree from CellSim 3D and the 3D view of the final cells are displayed. Colors has been assigned to the substances to make the 3D veiw match the input tree.

The tree has five final nodes and four internal nodes. 13 genes are used to code this tree. The complete DNA code is presented below:

# Start vector START: a0

# Internal node at level 1 g1(a1|a0)

g2(divide|a0)

# Internal nodes at level 2 g3(a2|a1,A)

g4(a3|a1,B) g5(divide|a1)

# Transition to final state g6(a4|a2,A)

g7(a5|a2,B) g8(a6|a3,B) g9(a2|a3,A)

(46)

34 CHAPTER 5. RESULTS g10(divide|a3,A) # Final state g11(a4|a4) g12(a5|a5) g13(a6|a6)

5.2.3 Simulating the tree of C. elegans

Here the tree of the nematode C. elegans is simulated. Figure 5.6 shows the simulation after 7 cell divisions. The cells are colored to show cells from three different subtrees. No simplifications are made when coding the tree, which means that about 4000 genes are used.

To generate a tree with 1024 final unique nodes there are 1023 internal nodes. Each internal node generates a gene of type gx(ax|ax−1, A/B). Also 1023 genes

of type gdx(divide|ax) are needed. 1024 genes of type gx(f inalx|ax − 1, A/B)

generates the final nodes and 1024 nodes of type gf inalx(af inalx|af inalx) keeps the

final state of a cell. This sums up to 1023 + 1023 + 1024 + 1024 = 4094 genes.

(47)

Chapter 6

Conclusions and future work

In this chapter the results are discussed and possible future work on the simulator is proposed.

6.1 Discussion

CellSim 3D provides the basic functionality for physical simulation of dividing cells. It has tools for visualizing the results of the simulation in an easy way. The simulator performs well up to around 128 cells, which corresponds to the first 7 cell divisions. This is enough for some purposes, but increased scalability whould be even better. An investigation of what the computational bottleneck is would be a good starting point for improvements of the simulator.

6.2 Future work

Some more features still remain to be implemented into the simulator. That includes cell to cell connections, and cell signaling. Also running more simulations in the existing simulator would be interesting. In this section some specially interesting investigations are mentioned.

6.2.1 Improvements of the physcical simulator

Lowering the resolution of the mesh object when the number of cells increases and the size of each cell decreases would improve the performance a bit. However the implemented collision handling system is designed to work best with equally sized objects with a mesh resolution that is about that used in the current version.

The wire frame view is a good tool when studying where cells or groups of cells are placed inside an organism. A better visualization with real transparancy

(48)

36 CHAPTER 6. CONCLUSIONS AND FUTURE WORK

would make it even easier to see this kind of information.

6.2.2 Cell connections

Cell to cell adhering connections and policies for making and releasing connections would be an interresting topic which have not been tested. Programming different policies and evaluate the resulting visual appearence would be interesting.

6.2.3 Cell placement

CellSim 3D gives the possibility to study how cells are placed inside an organism. An investigation of how cell placement is related to cell position in the lineage tree could be done using CellSim 3D. Also different policies for placing the divi-sion plane and how that affects the placement of cells could be included in this investigation.

6.2.4 Use physical data in the model of cell control

One possibility is to use the physical properties from the physical simulator to make decisions about cell differentiation and division. For example the type of neighbouring cells could affect how a cell differentiates and cell size could be a factor to include when a cell decides if it is about to divide or not.

6.2.5 Cell communication

Cell communication can be implemented to make cells differentiate the same way as the majority of its neighbours. This mechanism can probably make it possible to code the same organism with less complex code in some cases.

(49)

Bibliography

[1] B. Alberts, D. Bray, J. Lewis, M. Raff, K. Roberts, and J. D. Watson. Molecular Biology of the Cell. Third edition, 1994.

[2] Eric B. Becker, Graham F. Carey, and J. Tinsley Oden. Finite elements an intro-duction. 1981.

[3] Marie-Paule Cani-Gascuel and Mathieu Desbrun. Animation of deformable models using implicit surfaces. IEEE Transactions on visualization and computer graphics, 3(1), 3 1997.

[4] Bruno Heidelberger, Matthias Teschner, Richard Keiser, Matthias Mller, and Markus Gross. Consistent penetration depth estimation for deformable collision response, 2004.

[5] Hodgkin J. Introduction to genetics and genomics, September 2005.

[6] Tomas Jakobsen. Advanced character physics. In Proceedings of Game Developers Conference, 2001.

[7] Roger Johansson. Simulering av mjuka kroppar for spel, 2006.

[8] William E. Lorensen and Harvey E. Cline. Marching cubes: A high resolution 3d surface construction algorithm. Computer Graphics, 21(4), 7 1987.

[9] M. Matyka. How to implement pressure soft body model, 2004.

[10] M. Matyka and M. Ollila. A pressure model for soft body simulation. In Proc. of Sigrad, UMEA, 2003.

[11] David M. Prescott. Cells. Jones and Bartlett publishers, 1988. [12] Stephen Sorkin. Efficient collision detection for flexible objects, 2000. [13] J. Spieth and D. Lawson. Overview of gene structure, January 2006.

[14] Matthias Teschner, Bruno Heidelberger, Matthias Mller, Danat Pomeranets, and Markus Gross. Optimized spatial hashing for collision detection of deformable ob-jects, 2003.

[15] Gino van den Bergen. Efficient collision detection of complex deformable models using aabb trees, 1998.

(50)

(51)

Appendices

(52)

(53)

Appendix A

Calculating the volume of a mesh

The volume of a manyfold mesh can be calculated using the faces of the mesh. Gauss theorem A.1 relates the volume and surface integrals.

Z S F · dA = Z V ∇ · Fdτ (A.1)

The theorem states that the surface integral of a vector field F times the unit normal equals the volume integral of the divergence of the same vector field. Using a vector field with constant divergence ∇ · F = c gives equation A.2 that gives the volume times a constant c.

Z V ∇ · Fdτ = Z V cdτ = c Z V dτ = cV (A.2)

If F = (x, y, z), equationA.3-A.6gives the divergence 3. This menas that 3V is calculated by equation A.2.

∇ · F = ∇ · (x, y, z) (A.3) = ( δ δx, δ δy, δ δz) · (x, y, z) (A.4) = (δx δx, δy δy, δy δz) (A.5) = 3 (A.6)

EquationA.7approximates the integral in equationA.2with a Riemann sum. 3V = Z V ∇ · Fdτ = Z S F ≈X iǫS F(fi) × A(fi) (A.7) 41

(54)

42 APPENDIX A. CALCULATING THE VOLUME OF A MESH

In the approximation the vectorfield can be evaluated at each triangle, F(fi),

giving the formula in equationA.8, where v1, v2and v3are the vertices of triangle

i.

3V =X

iǫS

(v1+ v2+ v3)fi

Physical Simulation and Visualization of Cells

Contents

List of Terms

and Abbrevations

Chapter 1

Introduction

1.1

Goal

1.2

Outline

Chapter 2

Cells

2.1

Function and structure of cells

2.2

Cell types

2.3

The cell cycle and cell division

2.4

Genes and gene expression

2.5

Cell communication

2.6

Caenorhabditis elegans

Chapter 3

Soft body simulation

3.1

Modeling soft bodies

3.2

Particle systems

3.3

Implicit modeling

3.4

Finite elements method

3.5

Differential equations

3.6

Collision handling

Chapter 4

Method

4.1

Overview of the simulation

4.2

Soft body simulation

4.3

Collision handling

4.4

Cell simulation

Chapter 5

Results

5.1

CellSim 3D

5.2

Simulations

Chapter 6

Conclusions and future work

6.1

Discussion

6.2

Future work

Bibliography

Appendices

Appendix A

Calculating the volume of a mesh