**3. Cell-based modeling approaches**

As described in Sections 1 and 2, prion diseases offer a particularly intriguing biological phenomenon for mathematical and computational analysis because such diseases cover many different systems and spatial scales. At the level of a population, prion diseases can be studied as a classical epidemic model where infections are spread among a susceptible population. Alternatively, prion diseases can be studied on a microscopic scale as a genetic disease whose phenotype is caused by prion aggregates that cause a gain of function mutation in certain genes. Current models of protein aggregation in yeast have successfully provided further insight into important mechanisms driving prion disease dynamics (i.e., conversion and fragmentation), but there is a need to develop models that consider the underlying microscopic processes of protein aggregation together with macroscopic properties of the environment in which they are taking place. This requires modeling frameworks that consider the impact of processes taking place on many different spatiotemporal scales. One such class of models that has been developed primarily for studying biological processes from many different scales is cell-based models.

In this section we review cell-based modeling approaches and their various applications. Cell-based models are mathematical models that represent individual cells as discrete entities and produce simulations to predict large-scale, collective behavior of a population or group of cells from the behavior and interactions of individual cells. The inputs to a cell-based model are experimentally observed cell behaviors including how individual cells respond to both intracellular and extracellular cues. Cell behaviors are encoded in a set of biologically relevant rules for cells to follow during simulations. The outputs of a cell-based model are measurable population-level and cell-level characteristics that follow nontrivially from cell–cell coordination and individual cell response to the changing microenvironment, including changes caused by the cells themselves.

Within this class of models a further distinction is made between lattice-based or cellular automaton (CA) models in which the particles live on the coordinates of a lattice and off-lattice methods that use real numbers to describe the coordinates of each particle. Biologically motivated CA models describe cell-cell and cellenvironment interactions through phenomenological local rules, which allow for efficient simulation of many different biological systems ranging from bacteria, excitable media, and chemotactic aggregation to chicken embryonic tissues and tumors. (For reviews see [70–77].) CA models are well suited for studying the dynamics of systems with a large number of cells, i.e., growing monolayers or tumors, because their computation time is efficient and they allow for systematic sensitivity analysis of different parameters [70, 73, 75, 76, 78, 79]. In addition, CA models provide an easy starting point from which to derive continuum equations for cell densities which are particularly helpful in studying the mechanisms driving pattern formation and the growth process of multicellular systems with billions of cells [73, 78, 80, 81]. However, CA models are rule-based, and model parameters do not represent a direct physical description of biophysical properties of cells [73, 78]. For this reason, we focus our discussion in this chapter on lattice-free modeling approaches.

In lattice-free modeling approaches, individual cells are represented using either a single particle or agent or a group of particles or agents that are free to move to any location in the computational domain [75, 76, 78, 82–86]. In simulations, each biological cell is modeled as a discrete entity and endowed with well-defined individual characteristics such as intracellular reaction kinetics or detailed biophysical properties. As stated by T. Newman, cell shape and cell response to local mechanical forces is one of the most challenging biological characteristics to build into a cell model [78]. For this reason, lattice-free modeling approaches were developed as a way to explicitly model cell shape as well as investigate how individual cells respond to local mechanical forces and interactions with neighboring cells.

In lattice-free modeling approaches, interactions between agents are described using forces or potential functions, and the position of an individual cell is determined by solving an equation of motion for each agent that belongs to that cell. Moreover, it is standard to make the following two simplifying approximations for the equation of motion of each agent: (1) motion may be described by considering each vertex to be embedded in a viscous medium that applies a drag force on it with mobility coefficient *η*, and (2) inertia is vanishing [83, 87–90]. This leads to firstorder dynamics with the evolution of the position **x***i* of agent *i* to be determined by

$$
\eta \frac{d\mathbf{x}\_i}{dt} = \mathbf{F}\_i \tag{10}
$$

Lattice-free modeling approaches vary in complexity from those that track the center of each cell as a single point [91–93] to those that represent individual cells as collections of membrane and cytoplasm nodes to allow for more biologically relevant, emerging cell shapes [89, 94, 95]. Although a main advantage of CA methods has been their computational efficiency, as computational power continues to increase, lattice-free modeling approaches offer equally powerful tools for realistically modeling biological cell behaviors such as shape change, polarized cell growth, division, differentiation, and intercellular signaling dynamics in great detail without adding computational complexity. Lattice-free models have been used to investigate the impact of individual cell behaviors and make predictions about mechanisms controlling many different biological and pathological processes such as determining shape of a tissue, size of a cell colony, or biological function of an organ [75, 76, 78, 82, 84, 86, 96, 97]. In this section, we review three main categories of lattice-free modeling approaches, namely, center-based models, vertex models, and subcellular element models, as well as provide examples of how lattice-free model-

*Determining cell connectivity in center-based models. (A) In the overlapping spheres approach, each cell has an intrinsic radius, Ri, and two cells interact if their corresponding spheres overlap, i.e., if they are within a specified fixed distance from each other. (B) Triangulation involves computing a triangulation of the simulation domain using cell centers as nodes. Edges of the resulting mesh determine which cells are connected. (C) The Voronoi tessellation determines a region for each cell defined as a set of points in the plane that are closer to that*

*Multi-Scale Mathematical Modeling of Prion Aggregate Dynamics and Phenotypes in Yeast…*

*DOI: http://dx.doi.org/10.5772/intechopen.88575*

ing approaches have been used to study various biological problems.

gulation, or Voronoi tessellation as described in **Figure 3**.

the cells, then the force on cell *i* due to cell *j* is written as

Center-based modeling (CBM) approaches, or cell center models, track the center of each cell as a single point and can be classified by two main components, a definition of cell connectivity (**Figure 3**) and a definition of how cells interact or a force law (**Figure 4**) [76, 78, 83, 91]. There are several commonly used methods for determining cell connectivity, including the overlapping spheres approach, trian-

In center-based modeling approaches, the cell-cell interaction force is usually written as a function of the location of each cell centers and almost always acts in the direction of the vector connecting the two interacting cells [76, 78, 83, 91]. Let *ri* be the position of the cell center of cell *i*, and let **F***ij* denote the force on cell *i* due to cell *j*. If we assume that the force **F***ij* acts in the direction of the vector connecting

*ri* � *rj*

k*ri* � *rj*k (11)

**F***ij* ¼ *Fij*

**3.1 Center-based models**

**75**

**Figure 3.**

*cell center than any other.*

where **F***i*ð Þ*t* denotes the total force acting on vertex *i* at time *t*. The mobility coefficient *η* determines the timescale over which mechanical relaxation occurs and is often referred to as the damping coefficient.

*Multi-Scale Mathematical Modeling of Prion Aggregate Dynamics and Phenotypes in Yeast… DOI: http://dx.doi.org/10.5772/intechopen.88575*

#### **Figure 3.**

In this section we review cell-based modeling approaches and their various applications. Cell-based models are mathematical models that represent individual cells as discrete entities and produce simulations to predict large-scale, collective behavior of a population or group of cells from the behavior and interactions of individual cells. The inputs to a cell-based model are experimentally observed cell behaviors including how individual cells respond to both intracellular and extracellular cues. Cell behaviors are encoded in a set of biologically relevant rules for cells to follow during simulations. The outputs of a cell-based model are measurable population-level and cell-level characteristics that follow nontrivially from cell–cell coordination and individual cell response to the changing microenvironment,

Within this class of models a further distinction is made between lattice-based or cellular automaton (CA) models in which the particles live on the coordinates of a lattice and off-lattice methods that use real numbers to describe the coordinates of

environment interactions through phenomenological local rules, which allow for efficient simulation of many different biological systems ranging from bacteria, excitable media, and chemotactic aggregation to chicken embryonic tissues and tumors. (For reviews see [70–77].) CA models are well suited for studying the dynamics of systems with a large number of cells, i.e., growing monolayers or tumors, because their computation time is efficient and they allow for systematic sensitivity analysis of different parameters [70, 73, 75, 76, 78, 79]. In addition, CA models provide an easy starting point from which to derive continuum equations for cell densities which are particularly helpful in studying the mechanisms driving pattern formation and the growth process of multicellular systems with billions of cells [73, 78, 80, 81]. However, CA models are rule-based, and model parameters do not represent a direct physical description of biophysical properties of cells [73, 78]. For this reason, we

each particle. Biologically motivated CA models describe cell-cell and cell-

focus our discussion in this chapter on lattice-free modeling approaches.

to local mechanical forces and interactions with neighboring cells.

*η d***x***i*

is often referred to as the damping coefficient.

**74**

where **F***i*ð Þ*t* denotes the total force acting on vertex *i* at time *t*. The mobility coefficient *η* determines the timescale over which mechanical relaxation occurs and

*dt* ¼ **F***i* (10)

In lattice-free modeling approaches, individual cells are represented using either a single particle or agent or a group of particles or agents that are free to move to any location in the computational domain [75, 76, 78, 82–86]. In simulations, each biological cell is modeled as a discrete entity and endowed with well-defined individual characteristics such as intracellular reaction kinetics or detailed biophysical properties. As stated by T. Newman, cell shape and cell response to local mechanical forces is one of the most challenging biological characteristics to build into a cell model [78]. For this reason, lattice-free modeling approaches were developed as a way to explicitly model cell shape as well as investigate how individual cells respond

In lattice-free modeling approaches, interactions between agents are described using forces or potential functions, and the position of an individual cell is determined by solving an equation of motion for each agent that belongs to that cell. Moreover, it is standard to make the following two simplifying approximations for the equation of motion of each agent: (1) motion may be described by considering each vertex to be embedded in a viscous medium that applies a drag force on it with mobility coefficient *η*, and (2) inertia is vanishing [83, 87–90]. This leads to firstorder dynamics with the evolution of the position **x***i* of agent *i* to be determined by

including changes caused by the cells themselves.

*Apolipoproteins,Triglycerides and Cholesterol*

*Determining cell connectivity in center-based models. (A) In the overlapping spheres approach, each cell has an intrinsic radius, Ri, and two cells interact if their corresponding spheres overlap, i.e., if they are within a specified fixed distance from each other. (B) Triangulation involves computing a triangulation of the simulation domain using cell centers as nodes. Edges of the resulting mesh determine which cells are connected. (C) The Voronoi tessellation determines a region for each cell defined as a set of points in the plane that are closer to that cell center than any other.*

Lattice-free modeling approaches vary in complexity from those that track the center of each cell as a single point [91–93] to those that represent individual cells as collections of membrane and cytoplasm nodes to allow for more biologically relevant, emerging cell shapes [89, 94, 95]. Although a main advantage of CA methods has been their computational efficiency, as computational power continues to increase, lattice-free modeling approaches offer equally powerful tools for realistically modeling biological cell behaviors such as shape change, polarized cell growth, division, differentiation, and intercellular signaling dynamics in great detail without adding computational complexity. Lattice-free models have been used to investigate the impact of individual cell behaviors and make predictions about mechanisms controlling many different biological and pathological processes such as determining shape of a tissue, size of a cell colony, or biological function of an organ [75, 76, 78, 82, 84, 86, 96, 97]. In this section, we review three main categories of lattice-free modeling approaches, namely, center-based models, vertex models, and subcellular element models, as well as provide examples of how lattice-free modeling approaches have been used to study various biological problems.

#### **3.1 Center-based models**

Center-based modeling (CBM) approaches, or cell center models, track the center of each cell as a single point and can be classified by two main components, a definition of cell connectivity (**Figure 3**) and a definition of how cells interact or a force law (**Figure 4**) [76, 78, 83, 91]. There are several commonly used methods for determining cell connectivity, including the overlapping spheres approach, triangulation, or Voronoi tessellation as described in **Figure 3**.

In center-based modeling approaches, the cell-cell interaction force is usually written as a function of the location of each cell centers and almost always acts in the direction of the vector connecting the two interacting cells [76, 78, 83, 91]. Let *ri* be the position of the cell center of cell *i*, and let **F***ij* denote the force on cell *i* due to cell *j*. If we assume that the force **F***ij* acts in the direction of the vector connecting the cells, then the force on cell *i* due to cell *j* is written as

$$\mathbf{F}\_{ij} = F\_{ij} \frac{r\_i - r\_j}{||r\_i - r\_j||} \tag{11}$$

where *Fij* is the signed magnitude of **F***ij*. The total force on cell *i* is then

$$\mathbf{F}\_i = \sum\_j \mathbf{F}\_{ij} \tag{12}$$

cell bond breaking [76, 78, 83, 91]. For simplicity, we consider the linear law

*Multi-Scale Mathematical Modeling of Prion Aggregate Dynamics and Phenotypes in Yeast…*

attraction between distant cells, it is reasonable to define a cutoff range for

where *k*1 is a stiffness parameter. Since the linear law implies there is a large

*Fij* ¼ *k*1*δij* for *δij* ≥*δmin*

Some other choices for force laws include the Johnson-Kendall-Roberts (JKR) model [83, 91], the Hertz model [98, 99], and harmonic-like interaction [100–102]. The JKR model describes the interaction between two isotropic

homogeneous spheres that are strongly adhesive, i.e., as soon as the two cells come into contact range, they immediately form a contact area of finite size by active deformation of the cell membrane. However, a unique property of the JKR model is that when two cells are pulled apart, they continue to interact for distances greater than the standard contact range but less than a specified distance (i.e.,

the JKR model, the Hertz model approximates each cell as a homogeneous sphere, but in the Hertz model, hysteresis behavior is not included. In the Hertz model, the potential interaction is represented as the sum of a repulsive and an attractive

provide a simple model for the short-range cell-cell interaction due to cell-cell adhesion and elastic deformation that approximate two adhesively interacting cells by cuboidal objects with a non-deformable core, linked by linear springs. Each of these definitions for cell-cell interaction force offers their own set of advantages

Center-based models have been used to analyze multicellular processes in tumors [91, 99, 100], intestinal crypts and epithelial tissues [85, 103], cell migration in extracellular matrix [85, 100], and tissue regeneration, growth, and organization [91, 99, 100, 103, 104]. In addition, several center-based models have been developed for studying macroscopic properties of yeast colonies [92, 93, 105]. (See Section 4 for more details.) Although the mechanical information that can be extracted from these models is rather approximate, various biophysical aspects in these problems have explained important mechanisms for proper cell organization

Two-dimensional vertex models are a well-known class of lattice-free, cell-based models that provide a more detailed description of cell shape [84, 87, 107, 108]. In vertex models, each cell membrane is represented as a polygon composed of vertices and edges that are shared between adjacent cells (**Figure 5A**). A set of rules or an equation of motion defines how each vertex moves, and collective movement of vertices leads to changes in cell shape over time (**Figure 5A**). Most often, the equation of motion is a force or potential function based on the current

0 for *δij* , *δmin:*

*ij*). This phenomenon is referred to as hysteresis behavior. Similar to

*ijep*). On the other hand, harmonic-like force laws

*Fij* ¼ *k*1*δij* (16)

(17)

given by

*δmin* . *δij* . *δc*

force (i.e., *Fij* ¼ *Fa*

and disadvantages.

**3.3 Vertex models**

**77**

*ijttr* þ *Fr*

**3.2 Application of center-based models in biology**

and speed of growth in tissues [78, 97, 106].

interacting cells such that

*DOI: http://dx.doi.org/10.5772/intechopen.88575*

where the sum is over all cells *j* connected to cell *i*. As mentioned above, this force is usually taken to be balanced by a viscous drag as the cells move, so that the equation of motion for *ri*, the position of the cell center of cell *i*, is:

$$
\gamma \frac{dr\_i}{dt} = \mathbf{F}\_i \tag{13}
$$

where *γ* is the viscosity coefficient which could, for example, represent adhesion between a cell and the underlying substrate. In simulations there are several choices for the numerical method used to update the position of each cell center, but the most simple is forward Euler discretization. Thus, the position of a cell center at time *t* þ Δ*t* is given by

$$r\_i(t + \Delta t) = r\_i(t) + \frac{\Delta t}{\mathcal{Y}} \mathbf{F}\_i. \tag{14}$$

The cell–cell interaction force, **F***ij*, is usually a function of the overlap, *δij*, between two interacting cells which is given by

$$\delta\_{\vec{\eta}} = \mathcal{R}\_i + \mathcal{R}\_j - ||\mathbf{r}\_i - \mathbf{r}\_j||\tag{15}$$

where *Ri* and *Rj* are the radii of cell *i* and *j*, respectively, and *ri* and *rj* are the locations of the centers of cell *i* and cell *j*, respectively. Different definitions of the contact area between cells are possible to further extend the model [76, 78, 83, 91]. There exists a wide range of force laws that have been used in the literature, ranging from simple linear laws to more complex nonlinear models that can incorporate nonhomogeneous properties such as cell-cell adhesion, cell-substrate adhesion, and

#### **Figure 4.**

*Center-based model (CBM). The CBM approach updates cell positions based on interactions with connected cells. We show an interaction between two cells where each cell is considered as a sphere with a given radius (although other types of connectivity are possible (see* **Figure 3***)). The left hand side shows the position of node* ri *at time* t *with corresponding radius* Ri *and position of node* rj *at time* t *with corresponding radius* Rj*. The position of each node is used to calculate the overlap δij between cell* i *and cell* j *and the resulting linear force* F*ij to be applied to each node. The equation used to update the position of node* ri *is given as a function of the force* F*ij, time step* Δt*, and damping coefficient γ. The right hand side shows the new positions of node* ri *and* rj *at time* t þ Δt *after each position has been updated. (See Section 3.1 for more information).*

*Multi-Scale Mathematical Modeling of Prion Aggregate Dynamics and Phenotypes in Yeast… DOI: http://dx.doi.org/10.5772/intechopen.88575*

cell bond breaking [76, 78, 83, 91]. For simplicity, we consider the linear law given by

$$F\_{\vec{\imath}\vec{\jmath}} = k\_1 \delta\_{\vec{\imath}\vec{\jmath}} \tag{16}$$

where *k*1 is a stiffness parameter. Since the linear law implies there is a large attraction between distant cells, it is reasonable to define a cutoff range for interacting cells such that

$$F\_{\vec{\eta}} = \begin{cases} k\_1 \delta\_{\vec{\eta}} & \text{for} \quad \delta\_{\vec{\eta}} \ge \delta\_{\text{min}} \\ \mathbf{0} & \text{for} \quad \delta\_{\vec{\eta}} < \delta\_{\text{min}} \end{cases} \tag{17}$$

Some other choices for force laws include the Johnson-Kendall-Roberts (JKR) model [83, 91], the Hertz model [98, 99], and harmonic-like interaction [100–102]. The JKR model describes the interaction between two isotropic homogeneous spheres that are strongly adhesive, i.e., as soon as the two cells come into contact range, they immediately form a contact area of finite size by active deformation of the cell membrane. However, a unique property of the JKR model is that when two cells are pulled apart, they continue to interact for distances greater than the standard contact range but less than a specified distance (i.e., *δmin* . *δij* . *δc ij*). This phenomenon is referred to as hysteresis behavior. Similar to the JKR model, the Hertz model approximates each cell as a homogeneous sphere, but in the Hertz model, hysteresis behavior is not included. In the Hertz model, the potential interaction is represented as the sum of a repulsive and an attractive force (i.e., *Fij* ¼ *Fa ijttr* þ *Fr ijep*). On the other hand, harmonic-like force laws provide a simple model for the short-range cell-cell interaction due to cell-cell adhesion and elastic deformation that approximate two adhesively interacting cells by cuboidal objects with a non-deformable core, linked by linear springs. Each of these definitions for cell-cell interaction force offers their own set of advantages and disadvantages.

#### **3.2 Application of center-based models in biology**

Center-based models have been used to analyze multicellular processes in tumors [91, 99, 100], intestinal crypts and epithelial tissues [85, 103], cell migration in extracellular matrix [85, 100], and tissue regeneration, growth, and organization [91, 99, 100, 103, 104]. In addition, several center-based models have been developed for studying macroscopic properties of yeast colonies [92, 93, 105]. (See Section 4 for more details.) Although the mechanical information that can be extracted from these models is rather approximate, various biophysical aspects in these problems have explained important mechanisms for proper cell organization and speed of growth in tissues [78, 97, 106].

#### **3.3 Vertex models**

Two-dimensional vertex models are a well-known class of lattice-free, cell-based models that provide a more detailed description of cell shape [84, 87, 107, 108]. In vertex models, each cell membrane is represented as a polygon composed of vertices and edges that are shared between adjacent cells (**Figure 5A**). A set of rules or an equation of motion defines how each vertex moves, and collective movement of vertices leads to changes in cell shape over time (**Figure 5A**). Most often, the equation of motion is a force or potential function based on the current

where *Fij* is the signed magnitude of **F***ij*. The total force on cell *i* is then

equation of motion for *ri*, the position of the cell center of cell *i*, is:

time *t* þ Δ*t* is given by

**Figure 4.**

**76**

between two interacting cells which is given by

*Apolipoproteins,Triglycerides and Cholesterol*

*γ dri*

*ri*ð Þ¼ *t* þ Δ*t ri*ð Þþ*t*

The cell–cell interaction force, **F***ij*, is usually a function of the overlap, *δij*,

where *Ri* and *Rj* are the radii of cell *i* and *j*, respectively, and *ri* and *rj* are the locations of the centers of cell *i* and cell *j*, respectively. Different definitions of the contact area between cells are possible to further extend the model [76, 78, 83, 91]. There exists a wide range of force laws that have been used in the literature, ranging from simple linear laws to more complex nonlinear models that can incorporate nonhomogeneous properties such as cell-cell adhesion, cell-substrate adhesion, and

*Center-based model (CBM). The CBM approach updates cell positions based on interactions with connected cells. We show an interaction between two cells where each cell is considered as a sphere with a given radius (although other types of connectivity are possible (see* **Figure 3***)). The left hand side shows the position of node* ri *at time* t *with corresponding radius* Ri *and position of node* rj *at time* t *with corresponding radius* Rj*. The position of each node is used to calculate the overlap δij between cell* i *and cell* j *and the resulting linear force* F*ij to be applied to each node. The equation used to update the position of node* ri *is given as a function of the force* F*ij, time step* Δt*, and damping coefficient γ. The right hand side shows the new positions of node* ri *and* rj *at time*

t þ Δt *after each position has been updated. (See Section 3.1 for more information).*

**F***i* ¼ X *j*

where the sum is over all cells *j* connected to cell *i*. As mentioned above, this force is usually taken to be balanced by a viscous drag as the cells move, so that the

where *γ* is the viscosity coefficient which could, for example, represent adhesion between a cell and the underlying substrate. In simulations there are several choices for the numerical method used to update the position of each cell center, but the most simple is forward Euler discretization. Thus, the position of a cell center at

> Δ*t γ*

*δij* ¼ *Ri* þ *R j* �k*ri* � *rj*k (15)

**F***ij* (12)

*dt* ¼ **F***i* (13)

**F***i:* (14)

*E* ¼ X *N*

*DOI: http://dx.doi.org/10.5772/intechopen.88575*

be chosen in a biologically relevant manner.

ficient Γ*α*, and preferred cell area *A*ð Þ 0

*Fi* ¼ X *α*∈*Si*

*A*ð Þ 0

**Figure 5**.

*α*¼1

*Kα* 2

associated with it, the force *Fi* may be written explicitly as

*α*

� �∇*iAα* þ Γ*αPα*∇*iPα* � � �X

and so did not account for cell neighbor rearrangements [87]. However, later models were formulated to simulate top-down dynamics of epithelial cells in a plane which required modeling cell neighbor rearrangements [107]. Namely, to accurately describe planar epithelial cell dynamics, cells in the model must be allowed to form and break bonds by changing connectivity among vertices as described in

where now the first sum runs over cells associated with vertex *i* and the second

Many of the earliest vertex models simulated cross-sections of epithelial tissues

Vertex models were initially developed to study the packing of bubbles in foams [109]. Similar to foams, cells in epithelial tissues are tightly packed and mechanically coupled to their neighbors by adhesion molecules along their common interfaces [84, 87, 107, 108]. In addition, cells exert forces onto each other and their environment. Since vertex models are well suited for modeling tightly packed cell ensembles with negligible intercellular space, the vertex modeling framework was later adapted to study two-dimensional packing and rearrangement of apical cell surfaces in planar epithelia [84, 87, 107, 108]. In these studies, epithelial cells are described by a planar two-dimensional network of vertices that defines the apical cell surfaces as polygons with straight interfaces between neighboring cells. The underlying assumption of two-dimensional apical vertex models is that the main forces acting to deform the cells are generated along the apical cell surfaces or can be effectively absorbed in the apical representation. The work of Odell et al. [110] first demonstrated the use of vertex models to study how spatial patterning in either active forces or passive mechanical properties may lead to tissue deformation. Subsequent studies by several groups have examined a variety of alternative patterns of force and material properties that can give rise to similar tissue deforma-

An important behavior of biological cells that does not affect the dynamics of foams but plays a large role in the structure of biological tissues is the ability of cells to grow, divide, and die. Model components representing cell division and death

require the modification of tissue connectivity which can add significant

�*Kα Aα* � *A*ð Þ 0

sum runs over vertices sharing an edge with it.

**3.4 Application of vertex models in biology**

tions [84, 87, 88, 107, 108, 111].

**79**

*Aα* � *A*ð Þ 0 *α* � �2

� �

*Multi-Scale Mathematical Modeling of Prion Aggregate Dynamics and Phenotypes in Yeast…*

þ Γ*α* 2 *P*2 *α*

where the area elastic modulus *Kα*, bond tension parameter Λ*ij*, perimeter coef-

In vertex models that use this energy function, the choice of preferred cell area

*α* depends on how cell growth is incorporated into the model. For a given energy function *U*, the force on each vertex is determined by the negative derivative of the energy with respect to the coordinates of that vertex *Fi* ¼�∇*iU*, where ∇*i* denotes the gradient operator evaluated at **x***i*. Computing the gradient of Eq. (18) and exploiting the fact that the movement of vertex *i* affects only the energy of the cells

þX *i*, *j*

*α* are all model parameters whose values must

*j*

Λ*ijlij* (18)

Λ*ij*∇*i*ł*ij* (19)

#### **Figure 5.**

*Vertex model. Vertex models depict each cell membrane as a polygon of vertices and edges that are shared between adjacent cells. At each time step, an equation of motion determines how each vertex moves, and collective movement of vertices over time results in changes in cell shape. (A) Example forces acting on vertex* i *due to interactions with shared vertices in cell α according to Eq. (18). (B–D) Cells in the model form and break bonds due to neighbor rearrangements. Here we describe three basic cell neighbor rearrangements accounted for by most vertex models [87, 107]. (B) T1 transitions occur when two vertices sharing a short edge (defined by a minimum distance much smaller than the average length of an edge) merge into a single vertex which is then reassigned to two new vertices, thus changing the local network topology. (C) T2 transitions occur when a given cell shrinks to an approximately zero area and is removed to represent apoptosis. (D) T3 transitions occur when the intersection of a vertex with an edge is avoided by replacing the approaching vertex with two new vertices. For more detail see Section 3.3.*

configuration of vertices, i.e., location and connection between pairs of vertices, as well as other geometrical features such as area or perimeter of cells [84, 87, 107, 108]. Many vertex models also incorporate rules to govern changes in connections between vertices to allow for rearrangements in cell neighbor relationships [107]. The precise equations of motion used differ between models and may be adapted to suit each particular biological problem. Below we review some of the common forms of the equation of motion and typical rules for rearrangements of vertex connectivity.

In most cases, the equation of motion for vertices is deterministic since a reasonable assumption for tightly packed cell aggregates is that stochastic motion of cells is mitigated by strong interactions between cells [84, 87, 107, 108]. One difference among vertex models in the literature lies in the definition of the force. The choice of form for the function **F***i* reflects which forces are thought to dominate cell mechanics for the system being studied. Some commonly modeled forces include tension or elastic forces, due to the combined action of a cell's actomyosin cortex [107] and adherens junctions [107], or pressure, due to hydrostatic pressure [87]. Note that the forces acting on each vertex may either be given explicitly, or else an energy function may be specified, whose gradient is assumed to exert a force on each vertex.

For a concrete example, we will consider an equation of motion built from an energy function given by Farhadifar et al. which has recently seen extensive use in modeling wing disk epithelia [87]. This energy function encodes constraints associated with the limited ability of cells to undergo elastic deformations, volume changes, and other movements due to adhesion to other cells. The energy function is defined by

*Multi-Scale Mathematical Modeling of Prion Aggregate Dynamics and Phenotypes in Yeast… DOI: http://dx.doi.org/10.5772/intechopen.88575*

$$E = \sum\_{a=1}^{N} \left( \frac{K\_a}{2} \left( A\_a - A\_a^{(0)} \right)^2 + \frac{\Gamma\_a}{2} P\_a^2 \right) + \sum\_{i,j} \Lambda\_{\vec{\eta}} l\_{\vec{\eta}} \tag{18}$$

where the area elastic modulus *Kα*, bond tension parameter Λ*ij*, perimeter coefficient Γ*α*, and preferred cell area *A*ð Þ 0 *α* are all model parameters whose values must be chosen in a biologically relevant manner.

In vertex models that use this energy function, the choice of preferred cell area *A*ð Þ 0 *α* depends on how cell growth is incorporated into the model. For a given energy function *U*, the force on each vertex is determined by the negative derivative of the energy with respect to the coordinates of that vertex *Fi* ¼�∇*iU*, where ∇*i* denotes the gradient operator evaluated at **x***i*. Computing the gradient of Eq. (18) and exploiting the fact that the movement of vertex *i* affects only the energy of the cells associated with it, the force *Fi* may be written explicitly as

$$F\_i = \sum\_{a \in \mathcal{S}\_i} \left( -K\_a \left( A\_a - A\_a^{(0)} \right) \nabla\_i A\_a + \Gamma\_a P\_a \nabla\_i P\_a \right) - \sum\_j \Lambda\_{\vec{\eta}} \nabla\_i \mathbf{l}\_{\vec{\eta}} \tag{19}$$

where now the first sum runs over cells associated with vertex *i* and the second sum runs over vertices sharing an edge with it.

Many of the earliest vertex models simulated cross-sections of epithelial tissues and so did not account for cell neighbor rearrangements [87]. However, later models were formulated to simulate top-down dynamics of epithelial cells in a plane which required modeling cell neighbor rearrangements [107]. Namely, to accurately describe planar epithelial cell dynamics, cells in the model must be allowed to form and break bonds by changing connectivity among vertices as described in **Figure 5**.

#### **3.4 Application of vertex models in biology**

Vertex models were initially developed to study the packing of bubbles in foams [109]. Similar to foams, cells in epithelial tissues are tightly packed and mechanically coupled to their neighbors by adhesion molecules along their common interfaces [84, 87, 107, 108]. In addition, cells exert forces onto each other and their environment. Since vertex models are well suited for modeling tightly packed cell ensembles with negligible intercellular space, the vertex modeling framework was later adapted to study two-dimensional packing and rearrangement of apical cell surfaces in planar epithelia [84, 87, 107, 108]. In these studies, epithelial cells are described by a planar two-dimensional network of vertices that defines the apical cell surfaces as polygons with straight interfaces between neighboring cells. The underlying assumption of two-dimensional apical vertex models is that the main forces acting to deform the cells are generated along the apical cell surfaces or can be effectively absorbed in the apical representation. The work of Odell et al. [110] first demonstrated the use of vertex models to study how spatial patterning in either active forces or passive mechanical properties may lead to tissue deformation. Subsequent studies by several groups have examined a variety of alternative patterns of force and material properties that can give rise to similar tissue deformations [84, 87, 88, 107, 108, 111].

An important behavior of biological cells that does not affect the dynamics of foams but plays a large role in the structure of biological tissues is the ability of cells to grow, divide, and die. Model components representing cell division and death require the modification of tissue connectivity which can add significant

configuration of vertices, i.e., location and connection between pairs of vertices, as

*Vertex model. Vertex models depict each cell membrane as a polygon of vertices and edges that are shared between adjacent cells. At each time step, an equation of motion determines how each vertex moves, and collective movement of vertices over time results in changes in cell shape. (A) Example forces acting on vertex* i *due to interactions with shared vertices in cell α according to Eq. (18). (B–D) Cells in the model form and break bonds due to neighbor rearrangements. Here we describe three basic cell neighbor rearrangements accounted for by most vertex models [87, 107]. (B) T1 transitions occur when two vertices sharing a short edge (defined by a minimum distance much smaller than the average length of an edge) merge into a single vertex which is then reassigned to two new vertices, thus changing the local network topology. (C) T2 transitions occur when a given cell shrinks to an approximately zero area and is removed to represent apoptosis. (D) T3 transitions occur when the intersection of a vertex with an edge is avoided by replacing the approaching vertex with two new vertices.*

107, 108]. Many vertex models also incorporate rules to govern changes in connections between vertices to allow for rearrangements in cell neighbor relationships [107]. The precise equations of motion used differ between models and may be adapted to suit each particular biological problem. Below we review some of the common forms of the equation of motion and typical rules for rearrangements of

In most cases, the equation of motion for vertices is deterministic since a reasonable assumption for tightly packed cell aggregates is that stochastic motion of cells is mitigated by strong interactions between cells [84, 87, 107, 108]. One difference among vertex models in the literature lies in the definition of the force. The choice of form for the function **F***i* reflects which forces are thought to dominate cell mechanics for the system being studied. Some commonly modeled forces include tension or elastic forces, due to the combined action of a cell's actomyosin cortex [107] and adherens junctions [107], or pressure, due to hydrostatic pressure [87]. Note that the forces acting on each vertex may either be given explicitly, or else an energy function may be specified, whose gradient is assumed to exert a force

For a concrete example, we will consider an equation of motion built from an energy function given by Farhadifar et al. which has recently seen extensive use in modeling wing disk epithelia [87]. This energy function encodes constraints associated with the limited ability of cells to undergo elastic deformations, volume changes, and other movements due to adhesion to other cells. The energy function

well as other geometrical features such as area or perimeter of cells [84, 87,

vertex connectivity.

*For more detail see Section 3.3.*

*Apolipoproteins,Triglycerides and Cholesterol*

**Figure 5.**

on each vertex.

is defined by

**78**

computational complexity, but these cell behaviors are essential mechanisms of many of the biological process being studied. Thus, one biological question that has been addressed using vertex models extended to include cell death and division is the control of packing geometries in an epithelial sheet [84, 87, 88, 107, 108, 111]. A highly cited example is the work of Farhadifar et al. [87] who performed a systematic analysis of the equilibrium cell packing geometries and their dependence on cell mechanical and proliferative parameters using the *Drosophila* wing disk as a model system. By comparing simulations with experimental results, the authors arrived at a set of parameter values for which their model accounts for vertex movements based on biologically relevant cell area variations, cell rearrangements due to laser ablation, and epithelial packing geometries seen in vivo. This work demonstrates how vertex models may be parametrized and tested using experimental data.

In many biological systems, patterns of mechanical stress and resulting macroscopic tissue shape changes may happen concurrently and can feed-back into each other. Vertex models can be easily modified to incorporate both mechanical and chemical feedbacks. For example, Odell et al. [110] included a simplified feedback in their earliest model by assuming that contraction of individual cells was activated by high levels of stretching due to tissue deformation. Although vertex models have been successful in testing many of the mechanisms that govern the macroscopic behavior of tightly packed epithelial tissues, this class of model generally ignores contributions from important components such as cell-matrix interactions [112], actomyosin contractility and other biophysical properties of the cell membrane [113], and active remodeling of cytoskeletal components. Thus, there have been several notable extensions of the vertex model to address these important factors [110, 114, 115]. Although vertex models provide a promising modeling framework for further investigation of different biological systems, in the next section, we review a different modeling framework that was developed to take into account more detailed representations of cell deformation and provide an extension to the types and complexity of mechanical stress patterning and feedback considered.

#### **3.5 Subcellular element model**

The subcellular element (SCE) model provides a framework for simulating lattice-free multicellular structures in which the shape of each cell dynamically emerges from interactions with the local environment due to model assumptions about the underlying mechanical properties of the system [82, 86, 89, 90]. In the model, each cell is composed of a large, and possibly varying, number of small nodes called subcellular elements (**Figure 6**). Each subcellular element of a cell is modeled as a single point at its center of mass, which changes position over time subject to three processes: (i) weak random fluctuations, (ii) elastic interaction with elements of the same cell, and (iii) elastic interaction with elements of other cells.

of neighboring segments of cytoskeleton) and medium-range attraction between elements of the same or different cells (modeling the adhesive forces between

*Subcellular element model (SCE). The SCE model uses sets of elements/nodes to represent individual cells. At each time step, an equation of motion is used to update the location of each node based on resulting forces from interactions between nodes due to cell growth, deformation, etc. (A) Illustration of a SCE model where membrane and cytoplasm of individual cells are represented by one set of homogeneous elements/nodes. Elements/nodes of the same cell are held together via intracellular interactions (solid lines) modeled using shortrange potential functions (*∇*Eintra). Elements/nodes of two neighboring cells interact via intercellular interactions (dashed line) modeled using potential functions that weakly bind adjacent elements/nodes (*∇*Einter). The force applied to node αi is calculated based on interactions with all other nodes βi* ∈*celli and adjacent nodes α j in neighboring cell, cellj. (B) SCE model where the membrane (filled nodes) and cytoplasm (clear nodes) of one cell are represented separately using two different sets of elements/nodes. The force applied to cytoplasm node αi is calculated based on intracellular interactions which include interactions with all other cytoplasm nodes in the same cell (EII) and interactions with all other membrane nodes of the same cells (EMI) via short-range potential functions. The force applied to membrane node βi is calculated based on both intracellular interactions and intercellular interactions. Intracellular interactions for membrane node β j include interactions with all cytoplasm nodes of the same cell (EMI) and interactions with adjacent membrane nodes of the same cell (EMMS), and intercellular interactions for membrane node β j include interactions with adjacent membrane nodes from neighboring cells. For more details about these interactions, refer to Section 3.5.*

*Multi-Scale Mathematical Modeling of Prion Aggregate Dynamics and Phenotypes in Yeast…*

*DOI: http://dx.doi.org/10.5772/intechopen.88575*

The potential functions described above are used in model equations to calculate the displacement of each element/node at each time step based on their interactions with neighboring elements/nodes resulting in the deformation of cells within the tissue (**Figure 6**). As mentioned previously, nodes are assumed to be in an overdamped regime so that inertia forces acting on the nodes are neglected [83, 87–90]. This leads to the following two equations of motion describing the

∇*E*intra �X

where *η* is the damping coefficient, *xi* is the position of node *i*, *N* is the total number of nodes, and ∇*E* is the potential function used to describe the interaction of node *i* with node *j*. This equation is discretized in time using the forward Euler method, and the position of node *xi* is incremented at discrete times as follows:

ð Þ¼ *t* þ Δ*t xαi*

*i*6¼ *j*

ðÞ�*t Fαi*

X *αi*

∇*E*inter

Δ*t η*

1

A (20)

(21)

segments of cytoskeleton) [116].

**Figure 6.**

movement of element/node *αi* in cell *i*:

where Δ*t* is the time step size.

**81**

*Fαi* ¼� X

*xαi*

0 @

*αi*6¼*βi*

In the SCE modeling approach, the membrane and cytoplasm of each cell can be represented together as one set of homogeneous elements/nodes or separately using two different sets of elements/nodes (**Figure 6**). Collective interactions between pairs of nodes from the same cell represent bulk cytoplasmic contents, and collective interactions between pairs of adjacent nodes from neighboring cells represent volume exclusion of cells as well as adhesive properties. In models that consider the membrane and cytoplasm separately, collective interactions between pairs of cytoplasm nodes and membrane nodes represent the cytoplasmic pressure of the cell applied to the cell membrane. Biomechanical and adhesive properties of cells are modeled through viscoelastic interactions between elements represented by phenomenological potential functions such as linear springs or Morse potential functions, which are used to simulate close-range repulsion (modeling volume exclusion *Multi-Scale Mathematical Modeling of Prion Aggregate Dynamics and Phenotypes in Yeast… DOI: http://dx.doi.org/10.5772/intechopen.88575*

#### **Figure 6.**

computational complexity, but these cell behaviors are essential mechanisms of many of the biological process being studied. Thus, one biological question that has been addressed using vertex models extended to include cell death and division is the control of packing geometries in an epithelial sheet [84, 87, 88, 107, 108, 111]. A highly cited example is the work of Farhadifar et al. [87] who performed a systematic analysis of the equilibrium cell packing geometries and their dependence on cell mechanical and proliferative parameters using the *Drosophila* wing disk as a model system. By comparing simulations with experimental results, the authors arrived at a set of parameter values for which their model accounts for vertex movements based on biologically relevant cell area variations, cell rearrangements due to laser ablation, and epithelial packing geometries seen in vivo. This work demonstrates how vertex models may be parametrized and tested using experimental data.

*Apolipoproteins,Triglycerides and Cholesterol*

In many biological systems, patterns of mechanical stress and resulting macroscopic tissue shape changes may happen concurrently and can feed-back into each other. Vertex models can be easily modified to incorporate both mechanical and chemical feedbacks. For example, Odell et al. [110] included a simplified feedback in their earliest model by assuming that contraction of individual cells was activated by high levels of stretching due to tissue deformation. Although vertex models have been successful in testing many of the mechanisms that govern the macroscopic behavior of tightly packed epithelial tissues, this class of model generally ignores contributions from important components such as cell-matrix interactions [112], actomyosin contractility and other biophysical properties of the cell membrane [113], and active remodeling of cytoskeletal components. Thus, there have been several notable extensions of the vertex model to address these important factors [110, 114, 115]. Although vertex models provide a promising modeling framework for further investigation of different biological systems, in the next section, we review a different modeling framework that was developed to take into account more detailed representations of cell deformation and provide an extension to the types and complexity of mechanical stress patterning and feedback considered.

The subcellular element (SCE) model provides a framework for simulating lattice-free multicellular structures in which the shape of each cell dynamically emerges from interactions with the local environment due to model assumptions about the underlying mechanical properties of the system [82, 86, 89, 90]. In the model, each cell is composed of a large, and possibly varying, number of small nodes called subcellular elements (**Figure 6**). Each subcellular element of a cell is modeled as a single point at its center of mass, which changes position over time subject to three processes: (i) weak random fluctuations, (ii) elastic interaction with elements of the same cell, and (iii) elastic interaction with elements of other cells. In the SCE modeling approach, the membrane and cytoplasm of each cell can be represented together as one set of homogeneous elements/nodes or separately using two different sets of elements/nodes (**Figure 6**). Collective interactions between pairs of nodes from the same cell represent bulk cytoplasmic contents, and collective interactions between pairs of adjacent nodes from neighboring cells represent volume exclusion of cells as well as adhesive properties. In models that consider the membrane and cytoplasm separately, collective interactions between pairs of cytoplasm nodes and membrane nodes represent the cytoplasmic pressure of the cell applied to the cell membrane. Biomechanical and adhesive properties of cells are modeled through viscoelastic interactions between elements represented by phenomenological potential functions such as linear springs or Morse potential functions, which are used to simulate close-range repulsion (modeling volume exclusion

**3.5 Subcellular element model**

**80**

*Subcellular element model (SCE). The SCE model uses sets of elements/nodes to represent individual cells. At each time step, an equation of motion is used to update the location of each node based on resulting forces from interactions between nodes due to cell growth, deformation, etc. (A) Illustration of a SCE model where membrane and cytoplasm of individual cells are represented by one set of homogeneous elements/nodes. Elements/nodes of the same cell are held together via intracellular interactions (solid lines) modeled using shortrange potential functions (*∇*Eintra). Elements/nodes of two neighboring cells interact via intercellular interactions (dashed line) modeled using potential functions that weakly bind adjacent elements/nodes (*∇*Einter). The force applied to node αi is calculated based on interactions with all other nodes βi* ∈*celli and adjacent nodes α j in neighboring cell, cellj. (B) SCE model where the membrane (filled nodes) and cytoplasm (clear nodes) of one cell are represented separately using two different sets of elements/nodes. The force applied to cytoplasm node αi is calculated based on intracellular interactions which include interactions with all other cytoplasm nodes in the same cell (EII) and interactions with all other membrane nodes of the same cells (EMI) via short-range potential functions. The force applied to membrane node βi is calculated based on both intracellular interactions and intercellular interactions. Intracellular interactions for membrane node β j include interactions with all cytoplasm nodes of the same cell (EMI) and interactions with adjacent membrane nodes of the same cell (EMMS), and intercellular interactions for membrane node β j include interactions with adjacent membrane nodes from neighboring cells. For more details about these interactions, refer to Section 3.5.*

of neighboring segments of cytoskeleton) and medium-range attraction between elements of the same or different cells (modeling the adhesive forces between segments of cytoskeleton) [116].

The potential functions described above are used in model equations to calculate the displacement of each element/node at each time step based on their interactions with neighboring elements/nodes resulting in the deformation of cells within the tissue (**Figure 6**). As mentioned previously, nodes are assumed to be in an overdamped regime so that inertia forces acting on the nodes are neglected [83, 87–90]. This leads to the following two equations of motion describing the movement of element/node *αi* in cell *i*:

$$F\_{a\_i} = -\left(\sum\_{a\_i \neq \beta\_i} \nabla E\_{\text{intra}} - \sum\_{i \neq j} \sum\_{a\_i} \nabla E\_{\text{inter}}\right) \tag{20}$$

where *η* is the damping coefficient, *xi* is the position of node *i*, *N* is the total number of nodes, and ∇*E* is the potential function used to describe the interaction of node *i* with node *j*. This equation is discretized in time using the forward Euler method, and the position of node *xi* is incremented at discrete times as follows:

$$
\varkappa\_{a\_i}(t + \Delta t) = \varkappa\_{a\_i}(t) - F\_{a\_i} \frac{\Delta t}{\eta} \tag{21}
$$

where Δ*t* is the time step size.

One of the important features of the SCE modeling approach is the ability to change parameters of the potential functions that are used to describe interactions between elements to calibrate model representations of biomechanical properties of a particular type of a cell directly using experimental data. More specifically, the SCE model can be used to perform in silico bulk rheology experiments on a single cell in order to scale the parameters such that the passive biomechanical properties of each cell are independent of the number of elements used to represent each cell [117]. As a result, SCE simulation output captures the underlying biomechanical properties of the real biological system being studied and remains relevant regardless of the choice of the number of elements used in the model.

mechanical model of individual platelet behavior that allowed the authors to examine the relationship between platelet stiffness and movement in the fluid. This facilitated the generation of new hypothesis about mechanisms platelets use to adhere to injured sites on the blood vessel walls. In the context of developmental biology, the primary application of the SCE to date has been a computational study of primitive streak formation by Sandersius et al. [86]. In addition, the SCE model has been applied to model development in both animal [95] and plant [94] systems.

*Multi-Scale Mathematical Modeling of Prion Aggregate Dynamics and Phenotypes in Yeast…*

In this chapter, we described the medical and scientific importance of studying misfolded protein diseases as well as provided a broad description of several classes of mathematical models that can be used to further investigate the underlying mechanisms governing protein aggregation and propagation in a multicellular system. In Section 2, we gave an overview of different modeling frameworks that have been developed and validated for studying protein aggregation in both yeast and mammalian systems. In Section 3, we described in detail various cell-based modeling frameworks that have been successfully used to gain insight about the impact of individual cell behaviors on macroscale properties of tumors, developmental tissues, and yeast colonies. However, our main goal in this chapter was to familiarize the reader with each class of model in order to facilitate further discussion about how these two types of models could be combined to develop a more complete representation of prion disease dynamics within an actively growing and dividing

While current models of protein aggregation in yeast have successfully provided further insight into important mechanisms driving prion disease dynamics (i.e., conversion and fragmentation), they have been unable to recreate a number of important physiological characteristics including variable phenotype induction rates that result in sectored phenotypes among a single yeast colony (see Section 1 for details). One reason for this may be that a large number of models representing intracellular dynamics of protein misfolding diseases were developed for studying protein aggregation dynamics in isolation and these models disregard the contribution of individual cell behaviors within the growing yeast colony as a possible mechanism governing prion disease dynamics. A major open question in prion biology is to understand how prion aggregates spread between cells within a whole colony or tissue. Experimental observations such as sectoring provide compelling data that transmission mechanisms other than what is addressed by current aggregation-only-models must play a role in the presence and persistence of prion disease phenotypes in yeast colonies, i.e., processes such as conversion, fragmentation, nucleation, and even enzyme-mediated fragmentation alone cannot entirely explain the spread of prion disease throughout a yeast colony. Thus, in order to test hypotheses about the impact of individual cell behaviors on the spread of misfolded

proteins, it is necessary to develop a novel modeling framework.

Developing a modeling framework for investigating prion disease dynamics within an entire yeast colony is challenging because it requires capturing the physical processes on an individual cell level that determine colony growth (i.e., budding and variable cell cycle length) as well as capturing the interplay of individual cell processes with protein aggregation dynamics (i.e., asymmetric protein distribution at the time of division, persistence of diseased phenotype/aggregate that was given to daughter while it grows to begin a new cell cycle, lineage-dependent protein propagation). Several models have already been developed to investigate physical mechanisms controlling patterns of yeast colony growth such as cell division

**4. Conclusion**

*DOI: http://dx.doi.org/10.5772/intechopen.88575*

yeast colony.

**83**

As indicated in Fletcher et al. [84], computational experiments follow a creepstress protocol in which a constant extensile force is applied to the end of an SCE cell whose opposite end is fixed. Before the extensile force is released, the strain is measured as the extension of the cell in the direction of the force relative to its initial linear size. In silico estimates of the viscoelastic properties of cells modeled using the SCE approach have been shown in many biological applications to agree with in vitro rheology measurements [117, 118]. This indicates that the simple phenomenological dynamics of the SCE modeling approach are enough to capture low to intermediate responses of cytoskeletal networks over short timescales (10s) [118]. Over longer timescales (100 s), cells respond actively to external stresses by undergoing cytoskeletal remodeling, and this phenomenon can be incorporated into the SCE modeling approach by inserting and removing subcellular elements of a cell in regions under high or low stress [86].

The generalized Morse potential functions implemented in SCE modeling approach are commonly used in physics and chemistry to model intermolecular interactions [119] and in biology to represent volume exclusion of neighboring regions of the cytoskeleton [94, 95, 120–125]. While it is difficult to associate specific potential functions directly with specific cytoskeletal components of cells, computational studies of bulk properties at the tissue level have suggested that the precise functional form of the potential used in modeling has a small impact on overall system dynamics [97, 117].

#### **3.6 Application of SCE models in biology**

The subcellular element (SCE) modeling approach has been successful in modeling mechanical properties of individual cells as well as their components and determining individual cell impact on the emerging properties of growing multicellular tissue as well as describing cellular interactions with mediums such as the extracellular matrix and fluids [76, 84, 86, 89, 90, 95, 96, 111, 117, 120, 122, 124– 127]. The general modeling framework was initially developed by Newman et al. [89] for simulating the detailed dynamics of cell shapes as an emergent response to mechanical stimuli. Recent applications of the SCE modeling approach show that it is flexible enough to model additional diverse biological processes such as intracellular signaling [122], cell differentiation [94], and motion of cells in fluid [124].

To date, the SCE has not been widely used to study biological processes outside the area of epithelial morphogenesis. Christley et al. [122] developed a model of epidermal growth on a basal membrane that incorporates cell growth through the addition on new elements and division by redistributing a cell's group of elements between two new daughter cells. This mode was a novel implementation because it was coupled to a subcellular gene network representing intercellular Notch signaling. The SCE modeling framework has also been coupled to a fluid flow model to simulate the attachment of platelets to blood vessel walls [124]. Using the SCE framework for modeling individual platelets in simulations provides a detailed

*Multi-Scale Mathematical Modeling of Prion Aggregate Dynamics and Phenotypes in Yeast… DOI: http://dx.doi.org/10.5772/intechopen.88575*

mechanical model of individual platelet behavior that allowed the authors to examine the relationship between platelet stiffness and movement in the fluid. This facilitated the generation of new hypothesis about mechanisms platelets use to adhere to injured sites on the blood vessel walls. In the context of developmental biology, the primary application of the SCE to date has been a computational study of primitive streak formation by Sandersius et al. [86]. In addition, the SCE model has been applied to model development in both animal [95] and plant [94] systems.


**2. Materials and Methods**

The initial input database of 150 synthetic chemicals included JWH synthetic cannabinoids, others synthetic cannabinoids and other substances of abuse. These designer drugs were divided into three classes referred to as "Class 1—JWH", "Class 2—non-JWH Cannabinoids" and "Class 3—Others". The group of positives contained 50 JWH synthetic cannabinoids, while the group of negatives included 100 compounds, i.e., 50 non-JWH cannabinoids and 50 other substances of forensic interest [5].

We used the quantitative structure–activity relationship (QSAR) method to estimate and predict the pharmacokinetics, drug-likeness and medicinal chemistry friendliness of each input compound by calculating a number of 300 molecular descriptors in terms of their physical and chemical properties, as well as 50 indices characterizing their chemical

**Citation:** Burlacu, C.M.; Burlacu, A.C.; Praisler, M. Sensitivity Analysis of Artificial Neural Networks Identifying JWH Synthetic Cannabinoids Built with Alternative Training Strategies and Methods. *Inventions* **2022**, *7*, 82. https:// doi.org/10.3390/inventions7030082

Academic Editor: Anastasios Doulamis

Received: 13 July 2022 Accepted: 29 August 2022 Published: 13 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

109

absorption, distribution, metabolism, excretion and toxicity activity (ADMET). The descriptors were selected from three blocks, i.e., topological, 3D-MoRSE (molecule representation of structure based on electron diffraction) and toxicity [6].

Only the first 150 most relevant descriptors were selected and used for the final computational and modelling stage. Hence, the input database was a matrix consisting of 150 samples × 150 variables. The shape, feature and target types of the data set, including the list of the computed and tested input molecular descriptors was presented in a previously published article [7].

The data set was divided into three subsets of samples: training, selection and testing. Hence, we used 90 training samples (60%), 30 selection samples (20%) and 30 testing samples (20%). To discover redundancies between the input variables, we used a correlation matrix, which represents a numerical value between −1 and 1 that expresses the strength of the relationship between two variables [8]. The types of layers the most frequently used in our classification model were the perceptron layer, the probabilistic layer and the scaling and bounding layers. The objective of the selection was to find the best-performing network architecture in terms of system sensitivity.

To avoid underfitting and overfitting, the neuron selection algorithm responsible for finding the optimal number of neurons in the networks was the Growing neurons algorithm [9]. The *Neural Designer* data science and machine learning platform was used to generate the mathematical expression represented by ANNs in order to export and incorporate them into the programming language *Python 3.10* in the so-called production mode.

Our general training strategy consisted of two different concepts, i.e., the loss index and the optimization algorithms. The error was the essential term in the loss expression. The most important errors that we estimated were the sum squared error, the mean squared error, the root mean squared error, the normalized squared error and the Minkowski error. We used the L1 and L2 regularization methods, which involve the sum of the absolute values of all parameters and the square sum of all the parameters in the ANN. The loss index was measured on the data set and could be represented as a hyper-surface with the parameters as coordinates (see Figure 1) [10].

**Figure 1.** Loss index viewed as a hyper-surface with parameters as coordinates.

In order to train the ANN, we generated a sequence of parameter vectors so that the loss index was reduced at each iteration of the algorithm.

#### **3. Results**

Five different optimization algorithms were used and compared, each with a variety of different calculation and storage requirements: gradient descent [11], conjugate gradient, quasi-Newton method, Levenberg-Marquardt algorithm [12] and adaptative linear momentum [13].

In order to scale the inputs, we calculated the following parameters: the minimum, the maximum, the mean and the standard deviation (see Table 1). The ANN architecture is presented in Figure 2 for version 1. The architectures of the following versions (2, 3 and 4) were also with one hidden layer perceptrons and had the same input and output layers as version 1. On the other hand, their hidden layers contained three (version 2), one (version 3) and six (version 4) nodes respectively.

**Table 1.** Values used for scaling the inputs for all ANN versions.

**Figure 2.** Network architecture, version 1: scaling layer (yellow), perceptron layer (blue), probabilistic layer (red).

We used the adaptive moment estimation (version 1), the Levenberg–Marquardt (version 2), the gradient descent (version 3) and the conjugate gradient optimization algorithms, as well as the growing neuron selection (all versions) method with the L1 (versions 2, 3 and 4) and L2 (version 1) regularization methods.

#### **4. Discussion**

The confusion matrices, calculated for each architecture and 30 testing samples, are presented in Tables 2–5 and the error results are highlighted in Table 6. The activation functions used were the hyperbolic tangent (version 1), the rectified linear (versions 2, 3 and 4) and Softmax (all versions).

**Table 2.** Confusion matrix for the analyzed ANN—version 1: 29 (96.7%) tested compounds were correctly classified and 1 (3.3%) was misclassified.

**Table 3.** Confusion matrix for the analyzed ANN—version 2: 27 (90.0%) tested compounds were correctly classified and 3 (10.0%) were misclassified.

**Table 4.** Confusion matrix for the analyzed ANN—version 3: 26 (86.7%) tested compounds were correctly classified and 4 (13.3%) were misclassified.

**Table 5.** Confusion matrix for the analyzed ANN—version 4: 27 (90.0%) tested compounds were correctly classified and 3 (10.0%) were misclassified.

**Table 6.** Data errors for the analyzed ANNs.

In order to test and compare the performances of the analyzed ANNs, we used the weighted average derived from the confusion matrix, i.e., the accuracy, the recall and the F1 score (see Table 7) [14].

**Table 7.** Classification metrics for the target variables of each analyzed ANNs.

#### **5. Conclusions**

In terms of the system performance, the results obtained for the four ANNs designed to recognize the class identity of JWH synthetic cannabinoids lead to the following conclusions:

1. Accuracy [(true positives + true negatives)/total instances]:

In comparison with the accuracy (93.3%) obtained for the initial ANN model presented in a previous article, the ANN—amended version 1—generated a higher score (96.7%), while the other three ANNs generated a lower score (86.7% for the amended version 3 and 90.0% for the amended versions 2 and 4).

2. Sensitivity, or true positive rate (true positives/positive instances):

All the ANN architectures (the initial ANN and its four amended versions), had an exceptional sensitivity in detecting JWH synthetic cannabinoids (class I); all the tested JWHs were recognized as such without exception, the rate of true positives (TP) being 100% in all cases.

3. Specificity, or true negative rate (specificity = true negative/negative instances):

Compared to the specificity (90.0%) of the initial ANN model, one of the four new architectures recorded a higher score (95.0% for the amended version 1); the other three were characterized by a lower specificity (81.8% for amended version 3, 86.3% for amended version 2 and 86.3% for amended version 4).

4. Error rate [accuracy = (false positives + false negatives)/total instances]:

While the initial ANN model was characterized by an error rate of 6.6%, its amended version 1 recorded a better (lower) error rate (3.3%), while the other three ANNs recorded a higher error rate (10.0% for amended version 2, 3 and 4).

Regarding the goodness-of-fit analysis, the best results were recorded for the ANN amended version 1 (99.97%), followed by the amended version 2 (99.94%), amended version 3 (98.41%) and amended version 4 (99.04%). Hence, we may conclude that the best performing ANN architecture was the one that included the following elements:

The results indicate that very good classification rates were obtained although the data set was complex. We intend to continue this study by applying various metaheuristic algorithms on these data sets and compare the results.

**Author Contributions:** Conceptualization, C.M.B. and A.C.B.; methodology, M.P.; software, A.C.B.; validation, M.P.; formal analysis, C.M.B.; investigation, C.M.B.; resources, A.C.B.; writing—original draft preparation, A.C.B.; writing—review and editing, C.M.B.; supervision, M.P.; project administration, M.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** An appreciation of the "Wiley Online Library" and the data science platform and machine learning "Neural Designer", important and useful tools, used in the construction of system architectures presented in this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**