\documentclass[a4paper]{article}

\usepackage[utf8]{inputenc}

\usepackage[style=numeric]{biblatex}
\usepackage[acronym, nomain]{glossaries}
\usepackage{hyperref}
\usepackage{todonotes}

\addbibresource{references.bib}

\makeglossaries{}
\newacronym{cis}{CiS}{Cells in Silico}
\newacronym{cpu}{CPU}{Central Processing Unit}
\newacronym{cpm}{CPM}{Cellular Potts Model}
\newacronym{ecm}{ECM}{Extracellular Matrix}
\newacronym{fem}{FEM}{Finite Element Method}
\newacronym{gpu}{GPU}{Graphics Processing Unit}
\newacronym{lbm}{LBM}{Lattice Boltzmann Method}
\newacronym{mcs}{MCS}{Monte-Carlo Step}
\newacronym{mpi}{MPI}{Message Passing Interface}
\newacronym{nastja}{NAStJA}{Neoteric Autonomous Stencil code for Jolly Algorithms}
\newacronym{sls}{SLS}{Standard Linear Solid}

\begin{document}

\title{Exposé}
\author{Paul Brinkmeier}
\date{July 2023}
\maketitle

\section{Introduction}

Computational models of cell behavior can be useful to simulate and reiterate experiments.
In addition, they can show us how well our understanding models reality.
A popular approach is the \acrfull{cpm}, where each cell is modeled as a set of connected pixels or voxels on a two- or three-dimensional lattice.
To simulate biological processes involving thousands of cells, large lattices are needed.
Due to the local nature of the computations involved, the \acrshort{cpm} lends itself well to distributed programming.

In order to be true to \emph{in vivo}/\emph{in vitro} findings, such \emph{in silico} models must take into account a multitude of factors influencing cell behavior.
One such factor is the interaction with the \acrfull{ecm} that cells exist in.
In this work we will focus on the viscoelasticity of the collagen networks in the \acrshort{ecm}.
We explore models of viscoelasticity that, similar to the \acrshort{cpm} itself, employ local interactions to model global effects.
In this manner the simulations remain parallelizable.
Additionally, we investigate the performance of our model using different implementations on both \acrshort{cpu}s and \acrshort{gpu}s.

\section{Research}

\subsection{The \acrfull{cpm}}

The \acrshort{cpm}~\cite{graner1992} models cells as sets of sites on a square lattice which are usually connected.
Each lattice site is assigned the integer cell ID of the cell it belongs to.

The behavior of the cells is regulated by the Hamiltonian $H$, which represents the energy of a particular lattice.
The Hamiltonian contains at the least a term for the adhesion between cells on the lattice, but can be extended by other terms such as cell volume or alignment.

To advance the \acrshort{cpm}, a \acrfull{mcs} is performed:
The cell ID of a random lattice site is changed to the cell ID of one of its neighbors and the difference in energy $\Delta H$ is calculated.
The update is always accepted if the energy decreases.
If the energy does not decrease, the update is accepted probabilistically, where greater increases are less probable.
Repeated \acrshortpl{mcs} minimize $H$.

From an implementor's perspective, the \acrshort{cpm} has a great advantage over other approaches to cell simulation such as agent-based modeling:
Since updates happen on a square lattice and changes in energy can be calculated locally, it lends itself well to distributed programming.
\acrfull{cis}~\cite{berghoff2020} is a parallel implementation of the \acrshort{cpm} based on the \acrfull{nastja} framework~\cite{berghoff2018}.
\acrshort{nastja} offers an abstraction layer for implementing stencil codes on the \acrfull{mpi}, making it possible to leverage large-scale parallelism for \acrshort{cis}.
\acrshort{nastja} divides the domain into blocks, computing each stencil locally, then performing an halo exchanges of the border regions.

\subsection{The \acrfull{ecm}}

The \acrshort{ecm} is the part of a tissue that surrounds the cells.
It provides their physical and biochemical environment, thereby influencing cell behavior~\cite{frantz2010}.

While the \acrshort{ecm} consists of a variety of components, we focus on a single essential component:
Fibrous collagen networks and their viscoelasticity.
\acrshort{ecm} viscoelasticity has been established as an important factor in cell behavior~\cite{chaudhuri2020}.
For example, the \acrshort{ecm} confines cells and restricts processes such as migration, spreading, growth and mitosis.
These processes also affect the \acrshort{ecm} and can lead to permanent deformation.
In turn, this deformation can have an influence on cell behavior, resulting in self-reinforcing effects (see for example~\cite{mierke2021}).

\subsection{Models of the \acrshort{ecm} in the \acrshort{cpm}}

In this section we list current approaches to modeling the \acrshort{ecm} in \acrshort{cpm} simulations.
We focus on approaches that explicitly model the plasticity of \acrshort{ecm} collagens.

\paragraph{Static Cell}

A starting point is to model the \acrshort{ecm} as a static cell.
In this model, a cell ID is chosen to represent the solid parts of the \acrshort{ecm}.
Cell-matrix interactions are regulated by the Hamiltonian just like cell-cell interactions.
\acrshort{ecm} lattice sites do not copy their neighbors and can not be copied by their neighbors during a \acrshort{mcs}.
Instead, simulations using this approach usually allow cells to degrade adjacent matrix sites over time.
This approach is used for example in~\cite{bauer2007}, where the \acrshort{ecm} is initialized by randomly placing fiber bundles across the domain and~\cite{scianna2013}, which investigates cell behavior in \acrshortpl{ecm} with regular patterns.

\paragraph{Hybrid \acrshort{cpm}-\acrshort{fem}}

An approach using a \acrfull{fem} is presented in~\cite{vanoers2014} and expanded upon in~\cite{rens2017, rens2019}.
Each lattice site is assigned a local directional strain on the \acrshort{ecm}.
Cells exert traction forces on the \acrshort{ecm} used to calculate the lattice strains by a \acrshort{fem}.
The hamiltonian of the \acrshort{cpm} is modified such that cells respond to the strain.

\paragraph{Hybrid \acrshort{cpm} and Molecular Dynamics Methods}

Another approach is presented in~\cite{tsingos2022}.
This work simulates matrix fibers using a bead-and-chain model.
Similar to the previous approach, the \acrshort{ecm} model is coupled with the \acrshort{cpm}.
However, in this work, cells interact with the \acrshort{ecm} only through a sparse subset of lattice sites.

\subsection{Lattice Models of Viscoelastic Materials}

The strain response of the collagen networks in the \acrshort{ecm} is not fully elastic.
It exhibits both elastic (spring-like) and viscuous (damper-like) behavior.
The behavior of such viscoelastic materials is modeled by serial or parallel configurations of springs and dampers~\cite{sengul2021, mierke2021}.
The most common configurations for describing viscoelastic solids are

\begin{itemize}
	\item the Maxwell model, consisting of a spring and a damper in series,
	\item the Kelvin-Voigt model, consisting of a spring and a damper in parallel
	\item and the Zener or \acrfull{sls} model, which extends either the Maxwell or the Kelvin-Voigt model by another spring.
\end{itemize}

Depending on the specific viscoelastic characteristics that are to be predicted, a particular model can be chosen.
In order to align the viscoelastic \acrshort{ecm} model with the \acrshort{cpm}, we researched approaches that model viscoelastic materials on square lattices.

\paragraph{Discrete Particle Method}

A model for viscoelastic solids is presented in~\cite{obrien2008}.
This work extends the discrete particle method for elastic solids presented in~\cite{toomey2000}.
It is based on a two- or three-dimensional square lattice of particles.
Each particle is connected to all of its cardinal and diagonal neighbors.
The model for the force acting between two particles can be elastic or viscoelastic.
Various models are explored in~\cite{obrien2008, obrien2009, obrien2014, obrien2021}.

\paragraph{\acrfull{lbm}}

The \acrshort{lbm} is an established approach for modeling the dynamics of fluids~\cite{krueger2017}.
Also based on a square lattice, this model discretizes the particles moving at a particular lattice space into the cardinal and diagonal directions.
Research suggests that the \acrshort{lbm} can be used for modeling both solids~\cite{maquart2022} and viscoelastic fluids~\cite{malaspinas2010}.
Perhaps for this particular use case, a \acrshort{lbm} could be configured to model the \acrshort{ecm}.

\section{Contribution}

In this work we will explore lattice-based viscoelastic simulations of the \acrshort{ecm} in the \acrshort{cpm}.

\subsection{Method}

For the \acrshort{cpm} we use the distributed implementation \acrshort{cis}.
\acrshort{cis} is based on the \acrshort{nastja} framework implemented using the \acrshort{mpi}, which we will use to develop our model of the \acrshort{ecm}.
The model will likely be based on the viscoelastic discrete particle method, as preliminary experiments based on a simplified implementation show promising results.

In order to model cell-matrix interactions, we will develop a method that allows cells to influence the \acrshort{ecm} simulation.
To model matrix-cell interactions, we will expand the Hamiltonian of the \acrshort{cpm} to include a term dependent on the local configuration of the \acrshort{ecm}.
This should make it possible for our model to simulate the self-reinforcing interactions of cells and \acrshort{ecm}.

\subsection{Challenges}

Our preliminary experiments have produced some questions and likely challenges that our work will need to address.

\paragraph{Spatial Scale}

While we could simply use the same lattice for the \acrshort{ecm} model as for the \acrshort{cpm}, it is not clear that this will deliver the best results.
It could be useful to use a scaled lattice, e.g.\ where the lattice spacing of the \acrshort{ecm} model is twice as long.

\paragraph{Temporal Scale}

Compared to cells, the waves in a viscoelastic material move quickly.
It is likely that our model of the \acrshort{ecm} will have to go through multiple time steps between the \acrshortpl{mcs} of the \acrshort{cpm}.
In the context of \acrshort{nastja}, this means an increased number of halo exchanges between ranks per \acrshort{mcs}.
In order to reduce the number of halo exchanges, we could increase the width of the halo which allows the \acrshort{ecm} simulation to run for multiple time steps between halo exchanges.
As this approach necessarily leads to diminishing returns as the halo data gets bigger, an efficient configuration needs to be investigated.

\paragraph{Implementation Performance}

As \acrshort{cis} is designed to large and therefore compute-heavy simulations, it is worthwhile to measure the and optimize the compute needed by our implementation.
Since the discrete particle method is a dense approach, it should be possible to leverage common parallelization techniques such as vectorization and \acrshort{gpu} programming to improve performance.
In particular, it might prove useful to run the \acrshort{cpm} on \acrshortpl{cpu} and the \acrshort{ecm} model of \acrshortpl{gpu}.
We will experiment with these techniques and evaluate the possible improvements.

\newpage

\printglossary[type=\acronymtype, nogroupskip]

\printbibliography{}

\end{document}
