# Three facets of mathematical cancer biology research

Yue Wang  Department of Computational Medicine, University of California, Los
Angeles, California, United States of America

###### Abstract

Cancer, as the uncontrollable cell growth, is related to many branches of
biology. In this review, we will discuss three mathematical approaches for
studying cancer biology: population dynamics, gene regulation, and
developmental biology. If we understand all biochemical mechanisms of cancer
cells, we can directly calculate how the cancer cell population behaves.
Inversely, just from the cell count data, we can use population dynamics to
infer the mechanisms. Cancer cells emerge from certain genetic mutations,
which affect the expression of other genes through gene regulation. Therefore,
knowledge of gene regulation can help with cancer prevention and treatment.
Developmental biology studies acquisition and maintenance of normal cellular
function, which is inspiring to cancer biology in the opposite direction.
Besides, cancer cells implanted into an embryo can differentiate into normal
tissues, which provides a possible approach of curing cancer. This review
illustrates the role of mathematics in these three fields: what mathematical
models are used, what data analysis tools are applied, and what mathematical
theorems need to be proved. We hope that applied mathematicians and even pure
mathematicians can find meaningful mathematical problems related to cancer
biology.

##  1 Introduction

Cancer, as a group of more than 100 types of diseases, is one of the main
causes of human death. It is the dream of countless scientists to further
understand cancer and finally cure cancer. However, many experiments might
cause harm [36], and cannot be performed on human beings. Therefore, given
limited types of data, we have to develop indirect methods to study cancer
biology, especially with mathematical tools.

In this review, we will go through some mathematical biology studies in these
fields, and discuss how research in these fields improves our understanding of
cancer biology. This review is targeted to researchers who have backgrounds in
mathematics and want to contribute to cancer biology research. We will discuss
what mathematical models are used, what data analysis tools are applied, and
what mathematical theorems need to be proved.

Population dynamics studies how the population level changes along time under
different circumstances. The growth of tumor can be described by population
dynamics. In Section 2, we introduce some mathematical studies in cancer cell
population dynamics in normal environment and under treatment. These studies
illustrate what mathematical problems arise in population dynamics, and what
insights can be obtained from mathematical analyses of population data.
Specifically, there are still open problems arising from related stochastic
models.

Gene regulation describes how the expression of one gene affects the
expression of another gene. Cancer-related genetic mutations can deviate the
normal cell fate through gene regulation. Nevertheless, gene regulation is
difficult to determine directly by experiments. In Section 3, we introduce
some mathematical studies in the inference of gene regulation (both mutual
regulation and autoregulation) from certain types of data. These studies
represent a broad field of developing mathematical inference methods in gene
regulation and other subjects. Various problems in probability and discrete
mathematics need to be solved.

Developmental biology studies how organized normal tissues arise from a single
zygote. It is the opposite of cancer biology, which studies the emergence of
abnormal tissues. A developing embryo can invert the fate of cancer cells. In
Section 4, we introduce some mathematical studies in developmental biology,
which could provide insights of studying cancer biology from the opposite
direction. Specifically, we introduce an inference method for certain
experiments, which leads to a coloring problem in lattice.

There have been many papers and books discussing the role of mathematics in
cancer biology [37, 68, 1, 19, 2]. We hope that this review can provide a
little more insights and encourage applied mathematicians and even pure
mathematicians to work on interesting problems related to cancer biology and
contribute to the cure of cancer.

##  2 Population dynamics

The most significant feature of cancer cells is the uncontrollable growth.
Therefore, measuring and studying the growth of cancer cell population should
be a central topic in cancer research. Many studies build models for cancer
cells living in a free environment [28]. Some other studies consider the
effect of nutrients [25] and drugs [3].

Traditional population dynamics models are deterministic, based on ordinary
differential equations (ODEs) [78] or partial differential equations (PDEs)
[56]. Since cell growth is highly stochastic, various stochastic processes are
used to model cellular clonal growth and evolution: birth-death processes
[52], Poisson processes [7], Markov chains [21], random sums of birth-death
processes [16], and branching processes [28]. Other stochastic approaches,
such as random dynamical systems [74], lifted Markov chains [61], continuous
Markov processes [45, 73], or more complicated models based on the above
methods [79], also have the potential of modeling cell population growth.

###  2.1 Phenotypic equilibrium phenomena

Mathematical models based on different biological assumptions can produce
different predictions, which can be verified by experiments. For instance, in
the traditional hierarchical model [47, 29, 14], the transition from cancer
stem cells to non-stem cancer cells is irreversible. In recent studies, there
is evidence indicating that non-stem cancer cells can convert back to cancer
stem cells in certain cancer types [39, 11, 46, 72, 17]. In mathematical
models that consider such reversible transitions, there are two phenomena that
do not exist in the hierarchical model: overshoot and phenotypic equilibrium
[40, 13]. Overshoot means that if the initial population value is lower than
the steady state, the population might first increase above the steady state
and then return to it. Phenotypic equilibrium means that starting from any
proportions of different phenotypes, the final proportions always converge to
the same constants. For example, breast cancer has three cell phenotypes:
stem-like, basal, and luminal. Starting from any one phenotype, after
cultivating for some time, the final proportions are always stem-like $2\%$,
basal $97\%$, luminal $1\%$ [21]. Therefore, one can perform corresponding
experiments to detect the existence of such phenomena, and determine whether
the cancer cell type in the experiment has conversions from the non-stem state
to the stem state.

In the following, we will review different models for cancer cell population
dynamics, and corresponding explanations for the phenotypic equilibrium
phenomena. We consider $n$ phenotypes in the cancer cell population:
$Y_{1},Y_{2},\cdots,Y_{n}$. Define
$\vec{X}(t)=[X_{1}(t),X_{2}(t),\cdots,X_{n}(t)]$ to be the population at time
$t$, where $X_{i}(t)$ is the population of phenotype $Y_{i}$.
$P_{i}(t)=X_{i}(t)/[\sum_{i=1}^{n}X_{i}(t)]$ is the proportion of phenotype
$Y_{i}$, which is not defined if $\sum_{i=1}^{n}X_{i}(t)=0$. Define
$\vec{P}(t)=[P_{1}(t),P_{2}(t),\cdots,P_{n}(t)]$.

####  2.1.1 Markov chain model

The first explanation for the phenotypic equilibrium phenomena is based on a
Markov chain model [21]. In this model, $\vec{P}(t)$ satisfies

$\mathrm{d}\vec{P}(t)/\mathrm{d}t=\vec{P}(t)\mathbf{Q},$

where $\mathbf{Q}$ is the transition rate matrix, satisfying
$\vec{1}\mathbf{Q}^{\prime}=\vec{0}$. This is the evolution equation for a
Markov chain. From classical Markov chain theory, we have the following [41]:
If $\mathbf{Q}$ is irreducible, then starting from any initial probability
distribution $\vec{P}(0)$, $\vec{P}(t)$ always converges to the same vector
$\vec{v}$, which satisfies $\vec{v}\mathbf{Q}=\vec{v}$. This explains the
phenotypic equilibrium phenomena in the Markov chain model. However, this
model has been criticized as unrealistic [28], since cells can grow, and
different phenotypes might have different growth rates.

####  2.1.2 ODE model

A classic model for population dynamics is the linear ODE model. We assume
that the population (not the proportion) satisfies

$\mathrm{d}\vec{X}(t)/\mathrm{d}t=\vec{X}(t)\mathbf{A},$

where $\mathbf{A}=\\{a_{i,j}\\}$ for $i,j=1,\ldots,n$, and $a_{i,j}\geq 0$ for
$i\neq j$ [78]. To study the limiting behavior of $\vec{X}(t)$, we first need
a result in matrix theory [32]:

######  Theorem 1.

$\mathbf{A}$ has a real eigenvalue $\lambda_{1}$, such that for any eigenvalue
$\mu\neq\lambda_{1}$, $\Re(\mu)<\lambda_{1}$. $\lambda_{1}$ has a left
eigenvector $\vec{u}=(u_{1},\ldots,u_{n})$ that satisfies $u_{i}\geq 0$ and
$\sum_{i=1}^{n}u_{i}=1$.

Now we can present the convergence result in the ODE model [78]:

######  Theorem 2.

Assume that $\lambda_{1}$ is a simple eigenvalue of $\mathbf{A}$. Starting
from any initial value $\vec{X}(0)$ except for a zero-measure set, we have
$\lim_{t\to\infty}\vec{P}(t)\to\vec{u}$.

Therefore, the proportion $\vec{P}(t)$ will converge to the unique attracting
fixed point $\vec{u}$. This explains the phenotypic equilibrium phenomena in
the ODE model. Nevertheless, on single cell level, the cell growth is
essentially stochastic. The ODE model is just a mass action approximation.

####  2.1.3 Branching process model

A more realistic model for cell population dynamics is the multi-type
branching process model. In this model, different cells at the same time point
are independent. For each cell of the $Y_{i}$ type, the behavior is described
by the following:

$Y_{i}\stackrel{{\scriptstyle\alpha_{i}}}{{\to}}d_{i1}Y_{1}+d_{i2}Y_{2}+\cdots{}+d_{in}Y_{n}.$

Therefore, a $Y_{i}$ cell lives an exponential time with expectation
$1/\alpha_{i}$ and turns into $d_{i1}$ $Y_{1}$ cells, $d_{i2}$ $Y_{2}$ cells,
$\cdots$, $d_{in}$ $Y_{n}$ cells, where $d_{i1},d_{i2},\cdots{},d_{in}$ are
random variables in $\mathbb{N}$. The random coefficients $d_{i1},d_{i2},$
$\cdots{},d_{in}$ may not be independent, but they are independent with the
exponential waiting time. For example, a symmetric division $Y_{1}\to 2Y_{1}$
means $(d_{11},d_{12},\cdots,d_{1n})=(2,0,0,\cdots,0)$. A transition $Y_{1}\to
Y_{2}$ means $(d_{11},d_{12},\cdots,d_{1n})=(0,1,0,\cdots,0)$. A death
$Y_{1}\to\emptyset$ means $(d_{11},d_{12},\cdots,d_{1n})=(0,0,0,\cdots,0)$.
Define

$a_{i,i}=\alpha_{i}(\mathbb{E}d_{ii}-1)\geq-\alpha_{i},\ \ \
a_{i,j}=\alpha_{j}\mathbb{E}d_{ij}\text{ for }i\neq j.$

Then the matrix $\mathbf{A}=\\{a_{i,j}\\}$ corresponds to the coefficient
matrix in the ODE model. We have the same $\lambda_{1}$ and $\vec{u}$ as in
Theorem 1.

In the branching process model, due to randomness, we need to assume that no
phenotype $Y_{i}$ has $X_{i}(t)=0$ for any $t>T$. This condition is called
non-extinction. We can prove the convergence of $\vec{P}(t)$, namely a law of
large numbers, conditioned on non-extinction. In a monograph on branching
process in 1972, a law of large numbers is proven for certain $\mathbf{A}$
[4]:

######  Theorem 3.

Assume that $\lambda_{1}$ is positive, and it is a simple eigenvalue of
$\mathbf{A}$. Assume $\mathbf{A}$ is irreducible. Conditioned on non-
extinction, $\vec{P}(t)\to\vec{u}$ almost surely.

This law of large numbers can be used to explain the phenotypic equilibrium
phenomena. Nevertheless, the irreducibility of $\mathbf{A}$ means that any two
phenotypes can convert between each other (might through some other
phenotypes). This might not hold for some cancer types. The irreducibility
assumption is partially lifted in a 2004 paper [26] (only assume that the
phenotype with the highest growth rate can convert to any other phenotype),
and fully lifted in a 2017 paper [28]:

######  Theorem 4.

Assume that $\lambda_{1}$ is positive, and it is a simple eigenvalue of
$\mathbf{A}$. Conditioned on non-extinction, $\vec{P}(t)\to\vec{u}$ almost
surely.

This explains the phenotypic equilibrium phenomena in the branching process
model.

In the branching process model, the waiting time for cell division is
exponential, which is not realistic. To overcome this flaw, one can generalize
the law of large numbers to non-Markovian branching processes. The trick is to
add auxiliary phenotypes, so that the overall waiting time is close to the
real distribution. Nevertheless, this trick is not rigorous, and general cases
are still open. This illustrates that the need in cancer research might
inspire new mathematical results.

###  2.2 Inferring the existence of multiple phenotypes

The experiments for phenotypic equilibrium need to measure the cell states. In
fact, even without knowing the cell state, but just counting cell numbers, we
can infer the existence of multiple phenotypes. In some experiments, leukemia
cells were cultivated. Some populations started from 10 cells, while other
populations started from 1 cell. The cell number of each population was
measured everyday. By comparing cell population levels, we can find that after
reaching the same population level, the number of initial cells still affects
cell behavior. Specifically, when reaching cell area 5 (approximately 2500
cells), the growth rates of populations growing from 10 initial cells are
higher than those of populations from 1 initial cell. See Fig. 1 for
translated growth curves. This difference indicates the existence of multiple
phenotypes with different growth rates: starting from 10 initial cells, it is
very likely that at least one cell is of the fast type, which will dominate
the population; starting from 1 initial cell, the population might only have
the slow type, which corresponds to a lower growth rate [67]. See Fig. 2 for a
schematic illustration. Certainly, this explanation needs to be verified by
further experiments that measure the phenotypes, especially with gene
sequencing. Without such data, one needs to run a parameter scan to confirm
that the model is robust.

Figure 1: Growth curves of cell populations starting from the time of reaching
5 area units, with a logarithm scale for the $y$-axis. The $x$-axis is the
time from reaching 5 area units. Red (left) and blue (right) curves correspond
to 10 or 1 initial cell(s). Although starting from the same population level,
patterns are different for distinct initial cell numbers. The $N_{0}=1$-cell
group has higher diversity and lower average population level. This figure
corresponds to the contents of Subsection 2.2.  Figure 2: Schematic
illustration of the qualitative argument: Three cell types and growth patterns
(three colors) with different seeding numbers. One $N_{0}=10$-cell well will
have at least one fast type cell with high probability, which will dominate
the population. One $N_{0}=1$-cell well can only have one cell type, thus in
the microculture ensemble of replicate wells, three possible growth patterns
for wells can be observed. This figure corresponds to the contents of
Subsection 2.2.

###  2.3 Cancer cell under treatment

Various drugs can be used to kill cancer cells. Besides traditional chemical
drugs, aptamers [60] and T cells [15] can also be used. Therefore, the
population dynamics of cancer cells under treatment is also worth studying.

The behavior of cancer under treatment can be complicated. The tumor cell
regrowth after treatment is generally regarded as the result of a Darwinian
somatic evolution: Given sufficient genetic variability in a sufficiently
large initial cell population, it is very likely that the population contains
cells with genomic mutations that lead to drug-resistance. A single cell with
such a mutation will survive the treatment and clonally expand, thus driving
the tumor regrowth under treatment. If this simple Darwinian model holds, then
the recurrence time of cancer should be positively related to the drug dose.
The reason is that at a higher drug dose, the number of survived cells is
smaller, thus taking a longer time to recover. However, in some experiments,
it is observed that the treatment effect does not always become better when
the drug dose increases. The treatment effect gradually saturates. If the drug
dose is sufficiently high, the cancer cell population recovers even faster.
Due to such phenomena, we have to introduce more complicated mechanisms. One
possible approach is to assume that at higher drug doses, the killing rate
does not further increase, but the drug-induced transition rate from the
sensitive state to the resistant state still increases [3]. Therefore, when we
increase the drug dose, the treatment effect first increases, since the
killing rate increases. Then the killing rate is almost saturated, and the
treatment effect does not further increase. Finally, the killing rate stops at
its maximum, and the transition rate still increases, which decreases the
treatment effect.

Besides the phenomena that higher drug dose leads to worse treatment effects,
decreasing the drug dose might lead to overcompensation, which means that the
population recovers to a higher level than the initial population [70].
Therefore, the treatment of cancer is a highly complicated decision process.
It is necessary to introduce reinforcement learning (especially with
stochasticity) into this field [50, 65, 76, 66, 77].

##  3 Gene regulation

Coding genes are transcribed to mRNAs and translated to proteins. Genes
(through their corresponding proteins) can regulate the expressions of other
genes. The origin of cancer is certain genetic mutations (e.g., p53) [22].
Such mutations change the gene expression patterns through gene regulation
[49], and the cell landscape also transforms [35]. Therefore, the study of
gene regulation is important in understanding and curing cancer [5].

###  3.1 Framework for gene regulation inference

Genes and their regulatory relations form gene regulatory networks (GRNs).
Since the gene expression and regulation are mostly confined within living
cells, it is almost impossible to determine gene regulation with direct
biochemical methods. Nevertheless, there have been many methods to infer the
GRN structure, namely whether one gene regulates another, from experimental
data [6, 9, 23]. Such inference methods use models based on different
assumptions, and need different types of data. To unify different approaches
for this problem, there is a framework for data classification, which has four
major dimensions: (1) whether measuring gene expression levels, or the levels
of certain phenotypes that are determined by such genes; (2) whether the
measurement is for a single cell (stochastic) or for a large population of
cells (deterministic); (3) whether the measurement is for one time point or
multiple time points; (4) whether one can add interventions on certain genes.
If the measurement is on single-cell level and at multiple time points,
depending on whether the same cell can be measured multiple times, there is an
extra dimension of whether we have the joint probability distribution or just
marginal probability distributions. In total, there are 20 different data
types [63]. Different data types need different inference methods, and some
data types (especially non-interventional) are not useful in inferring gene
regulation. Due to the large amount of related data, machine learning is a
promising future direction for inferring gene regulatory relations [24]. See
Table 1 for the classification of data types and corresponding methods that
can be applied to infer gene regulations.

| One-Time | Time Series
---|---|---
| Non-
---
Intervention
Intervention |  | Non-
---
Intervention
Intervention
Gene Expression |  | Single-
---
Cell
| Scenario 1
---
MF+DAG:
partial
| Scenario 2
---
PB: full
DAG: partial
MF+DAG:
full
| Scenario 3 a/b
---
3a Joint:
UC: full
3b Marginal:
MF+DAG:
partial
| Scenario 4 a/b
---
4a Joint:
UC: full
4b Marginal:
LS: full
PB: full
DAG: partial
MF+DAG: full
Bulk |  | Scenario 5
---
No
| Scenario 6
---
PB: full
DAG: partial
| Scenario 7
---
No
| Scenario 8
---
LS: full
PB: full
DAG: partial
Phenotype |  | Single-
---
Cell
| Scenario 9
---
No
| Scenario 10
---
PB: partial
| Scenario 11 a/b
---
No
| Scenario 12 a/b
---
PB: partial
LS+DAG:
partial*
PB+LS+DAG:
partial*
Bulk |  | Scenario 13
---
No
| Scenario 14
---
PB: partial
| Scenario 15
---
No
| Scenario 16
---
PB: partial
LS+DAG:
partial*
PB+LS+DAG:
partial*
Table 1: GRN structure inference with different data types: under what
assumptions, what structures can be inferred. There are 16 scenarios
classified by the following dimensions of data types: Gene Expression vs.
Phenotype; Single-Cell vs. Bulk; One-Time vs. Time Series; Non-Interventional
vs. Interventional. In Scenarios 3/4/11/12, there is an extra dimension of
Joint vs. Marginal. There are different assumptions: PB: path blocking; DAG:
directed acyclic graph; MF: Markov and faithful; LS: linear system; UC:
unconditional. Full/partial/no means all/some/no GRN structures can be
inferred. For example, “MF+DAG: partial” means under MF assumption and DAG
assumption, GRN structure can be partially inferred. The asterisk * in
Scenarios 12/16 means for some identified regulatory relations, we cannot
determine whether they are activation or inhibition [63]. This table
corresponds to the contents of Subsection 3.1.

###  3.2 Inference methods for mutual regulation

We briefly introduce some inference methods for mutual regulation. Each method
may apply in multiple scenarios. Different methods need different assumptions.
DAG assumption means that the GRN is a directed graph that has no directed
cycles. MF assumption means that the joint probability distribution of gene
expressions provide exactly all and only the conditional independence
relations implied by the GRN. LS assumption means that the gene expression
satisfies a linear ODE system. PB assumption means that if all paths between
two genes in a GRN have been blocked by intervened genes, then adding
intervention on one gene cannot further affect the other gene.

####  3.2.1 Causal DAG method

This method needs the MF and DAG assumptions, and applies to Scenarios 1 and
3b [20]. In this case, the GRN is a causal DAG, and we can use the conditional
independence relations to determine the corresponding causal DAG. However,
some different causal DAGs (all sharing the same edges, but with different
directions) produce the same conditional independence relations. Therefore, we
can only partially determine the GRN structure.

####  3.2.2 Path blocking relation method in gene expression

This method needs the PB assumption, and applies to Scenarios 2, 4b, 6, 8
[63]. From such data types, we know whether some genes fully block all paths
between two genes. There is an edge $G_{1}\to G_{2}$ if and only if all other
genes cannot block the path from $G_{1}$ to $G_{2}$. All edges can be inferred
in this case.

####  3.2.3 Ancestor-descendant relation method

This method needs the DAG assumption, and applies to Scenarios 2, 4b, 6, 8
[63]. From such data types, we can use interventions to determine whether one
gene $G_{1}$ has a directed path to another gene $G_{2}$ in the GRN (ancestor-
descendant relations). We can determine an equivalent class that share the
same ancestor-descendant relations by the following theorem:

######  Theorem 5.

The following procedure describes how to determine certain edges with
ancestor-descendant relations. (1) If $G_{j}$ is not a descendant of $G_{i}$,
then we can determine that the edge $G_{i}\to G_{j}$ does not exist. (2) If
$G_{j}$ is a descendant of $G_{i}$, and $G_{i}$ has another descendant
$G_{k}$, which is an ancestor of $G_{j}$, then we cannot determine the
existence of the edge $G_{i}\to G_{j}$, since we can find two GRNs with the
same ancestor-descendant relations, one with $G_{i}\to G_{j}$, and one without
$G_{i}\to G_{j}$. (3) If $G_{j}$ is a descendant of $G_{i}$, and $G_{i}$ does
not have another descendant $G_{k}$, which is an ancestor of $G_{j}$, then we
can determine that the edge $G_{i}\to G_{j}$ exists.

Therefore, we can partially determine the GRN structure. About the number of
inferred edges by this method, we have a lower bound [63]:

######  Theorem 6.

If the GRN is a connected DAG with $n$ genes, then we can use ancestor-
descendant relations to identify at least $n-1$ edges.

####  3.2.4 Ancestor-descendant relation and causal DAG method

This method needs the MF and DAG assumptions, and applies to Scenarios 2 and
4b [63]. By the causal DAG method, we can determine all edges in the GRN, but
the directions of some edges are unknown. Using the ancestor-descendant
relation, we can determine such directions. Therefore, the GRN can be fully
determined.

####  3.2.5 Causal time series method

This method does not need any assumption, and applies to Scenarios 3a and 4a
[20]. This is basically the same as the causal DAG method, but the time series
data naturally determines the direction of a causal relation. Therefore, we
can use conditional independence relations to fully determine the GRN
structure.

####  3.2.6 Linear ODE method

This method needs the LS assumption, and applies to Scenarios 4b and 8 [43].
We have a linear system for the gene expression levels, each with the form

$\mathrm{d}X_{i}/\mathrm{d}t=\sum_{j=1}^{n}a_{ji}X_{j}.$

We can measure $X_{i}(t)$ at different time points under different
interventions, and calculate $\mathrm{d}X_{i}/\mathrm{d}t$. Then we obtain a
linear algebraic equation system for $a_{ji}$. We can solve this and determine
the GRN structure.

####  3.2.7 Path blocking relation method in phenotype

This method needs the PB assumption, and applies to Scenarios 10, 12a/b, 14,
16 [63]. From such data types, we know whether some genes fully block all
paths from a gene $G_{i}$ to a phenotype $P$. If a subset $\mathcal{S}$ of
$\\{G_{1},\ldots,G_{n}\\}\backslash\\{G_{i}\\}$ blocks $G_{i}$ to $P$, but any
proper subset of $\mathcal{S}$ cannot block $G_{i}$ to $P$, then $\mathcal{S}$
is called a minimal blocking set. If a blocking set $\mathcal{S}$ is not
minimal, then $\mathcal{S}$ contains a blocking subset that is minimal. Define
$\beta(G_{i})$ to be all minimal subsets that block $G_{i}$ to $P$. The
following theorem describes what edges can or cannot be inferred from the path
blocking relations [63].

######  Theorem 7.

The following procedure describes how to determine certain edges with path
blocking relations for the phenotype. (1) There is an edge $G_{i}\to P$ if and
only if $\beta(G_{i})=\emptyset$. (2) If there exists
$\mathcal{S}\in\beta(G_{i})$, so that $G_{j}\notin\mathcal{S}$, and
$\mathcal{S}$ cannot block $G_{j}$ to $P$, then there is no edge $G_{i}\to
G_{j}$. (3) If (2) is not satisfied, but there exists
$\mathcal{S}\in\beta(G_{i})$, so that $G_{j}\in\mathcal{S}$, then there is an
edge $G_{i}\to G_{j}$. (4) If $\beta(G_{i})=\emptyset$, or for any
$\mathcal{S}\in\beta(G_{i})$, we have $G_{j}\notin\mathcal{S}$, and
$\mathcal{S}$ blocks $G_{j}$ to $P$, then we cannot determine whether
$G_{i}\to G_{j}$ exists, since we can find two GRNs with the same path
blocking relations, one with $G_{i}\to G_{j}$, and one without $G_{i}\to
G_{j}$.

Using the above procedure, we can partially determine the GRN structure. About
the number of inferred edges by this method, we have a lower bound [63]:

######  Theorem 8.

If the GRN with $n$ genes has a directed path from each gene $G_{i}$ to the
phenotype $P$, then we can use path blocking relations for the phenotype to
identify at least $n$ edges.

###  3.3 Inference methods for autoregulation

Some genes can regulate (activate or inhibit) their own expressions, which is
called autoregulation. Certain types of data can be used to infer the
existence of gene autoregulation, such as time series data [48, 71, 51] or
interventional data for single-cell gene expression [27]. There also exists an
inference method that only requires one-time non-interventional single-cell
gene expression data [57]. Although different inference methods are based on
different types of models, and require different data types, the basic idea is
the same: build two models, one with autoregulation, and one without
autoregulation; then try to find different behaviors between these two models.
This idea can be generalized to other problems of detecting the existence of
certain mechanisms. If these two models are indistinguishable, it is not a
negative result either: one proves that this mechanism is not detectable in
this setting [10]! Here we introduce a simple Markov chain model for gene
expression and present an inference method for autoregulation.

Consider $m$ genes $G_{1},\ldots,G_{m}$ for a single cell. Denote their
expression levels by random variables $X_{1},\ldots,X_{m}$. Each state of this
Markov chain, $(X_{1}=n_{1},\ldots,X_{i}=n_{i},\ldots,X_{m}=n_{m})$, can be
abbreviated as $\boldsymbol{n}=(n_{1},\ldots,n_{i},\ldots,n_{m})$. For gene
$G_{i}$, the transition rate of $n_{i}-1\to n_{i}$ is $f_{i}(\boldsymbol{n})$,
and the transition rate of $n_{i}\to n_{i}-1$ is $g_{i}(\boldsymbol{n})n_{i}$.
Transitions with more than one step are not allowed. The master equation of
this process is

$\begin{split}\frac{\mathrm{d}\mathbb{P}(\boldsymbol{n})}{\mathrm{d}t}=&\sum_{i}\mathbb{P}(n_{1},\ldots,n_{i}+1,\ldots,n_{m})g_{i}(n_{1},\ldots,n_{i}+1,\ldots,n_{m})(n_{i}+1)\\\
+&\sum_{i}\mathbb{P}(n_{1},\ldots,n_{i}-1,\ldots,n_{m})f_{i}(\boldsymbol{n})\\\
-&\mathbb{P}(\boldsymbol{n})\sum_{i}[f_{i}(n_{1},\ldots,n_{i}+1,\ldots,n_{m})+g_{i}(\boldsymbol{n})n_{i}].\end{split}$

Define $\boldsymbol{n}_{\bar{i}}=(n_{1},\ldots,n_{i-1},n_{i+1},\ldots,n_{m})$.
Define $h_{i}(\boldsymbol{n})=f_{i}(\boldsymbol{n})/g_{i}(\boldsymbol{n})$ to
be the relative growth rate of gene $V_{i}$. Autoregulation means for some
fixed $\boldsymbol{n}_{\bar{i}}$, $h_{i}(\boldsymbol{n})$ is (locally)
increasing/decreasing with $n_{i}$, thus $f_{i}(\boldsymbol{n})$
increases/decreases and/or $g_{i}(\boldsymbol{n})$ decreases/increases with
$n_{i}$. In this model, we have the following result [57]:

######  Theorem 9.

In the above model, assume the GRN has no directed cycle that contains gene
$G_{k}$. Assume $g_{k}(\cdot)$ is a constant for all $\boldsymbol{n}$. If
$G_{k}$ has no autoregulation, meaning that $h_{k}(\cdot)$ and $f_{k}(\cdot)$
do not depend on $n_{k}$, then the expression level of $G_{k}$ satisfies
$\mathrm{var}(X_{k})/\mathbb{E}(X_{k})\geq 1$. Therefore, $G_{k}$ has
$\mathrm{var}(X_{k})/\mathbb{E}(X_{k})<1$ means $G_{k}$ has autoregulation.

Nevertheless, this method fails to determine whether autoregulation exists
when we have $\mathrm{var}(X_{k})/\mathbb{E}(X_{k})\geq 1$ for one gene. It is
hypothesized that the requirement for $g_{k}(\cdot)$ can be dropped:

######  Conjecture 1.

In the above model, assume the GRN has no directed cycle that contains gene
$G_{k}$. If $h_{k}(\cdot)$ and $f_{k}(\cdot)$ do not depend on $n_{k}$, then
the expression level of $G_{k}$ satisfies
$\mathrm{var}(X_{k})/\mathbb{E}(X_{k})\geq 1$.

###  3.4 Other topics in gene regulation

If we cannot add intervention, and cannot measure the same cell(s) multiple
times, then data from different time points are essentially the same, and we
cannot obtain information of the dynamics. In this case, sometimes it is
impossible to build the causal relationship in gene regulation (including
autoregulation) [62].

With the inferred knowledge of gene regulations, one can explain some
phenomena in cancer biology. For example, myeloproliferative neoplasm needs
two mutations, JAK2 and TET2. The order of the appearance of these two
mutations can affect clinical features and responses to therapy [42, 34, 33].
Possible explanations include threshold phenomena and hidden variables.

Gene regulation has many complicated mechanisms. Certain genes can translocate
in gene sequences, which are called transposons. Genes at different locations
can have different effects in gene regulation [38]. Therefore, it is also
important to design methods to determine transposons in gene sequences [55,
31, 18, 12, 75]. The general idea is to find the longest common subsequence,
whose complement is all the transposons. Determining the longest common
subsequence is a classic problem in computer science.

##  4 Developmental biology

Developmental biology studies how a zygote differentiates into normal tissues,
while cancer biology studies how normal cells become malfunctional. Seemingly
in opposite directions, these two subjects are interconnected. The
microenvironment of embryo can suppress the emergence of cancer [8]. In fact,
certain cancer cells transplanted to an embryo will differentiate into normal
tissues [44]. The development of an embryo is extremely complicated, so that a
cumbersome model is still an oversimplification [59].

###  4.1 Inference on transplantation experiments

The transplantation experiments of cancer cells and embryos indicate a
possible direction of cancer treatment. However, such transplantation
experiments are expensive and difficult to perform. Therefore, inference
methods that determine unknown experimental results from known results are
necessary [64]. The idea is inspired by Ising model in statistical physics.
Each transplantation experiment is determined by the donor tissue and the host
tissue. We can construct the similarities between tissues. Then
transplantation experiments can be compared by the similarities between their
donor tissues and host tissues. We can assume that similar experiments are
more likely to have similar results. Given some transplantation experiments,
where some results are known, and some results are unknown, we can take a
guess for those unknown results. Then we feed all results to a penalty
function. Assume the transplantation experiments are $X_{1},\ldots,X_{n}$;
their results are $x_{1},\ldots,x_{n}$; the similarity function between two
transplantation experiments is $F(X_{i},X_{j})$; the similarity function
between results is $f(x_{i},x_{j})$. Then the penalty function should be

$P(x_{1},\ldots,x_{n})=-\sum_{i=1}^{n}\sum_{j=1}^{n}F(X_{i},X_{j})f(x_{i},x_{j}).$

Therefore, for each result configuration $(x_{1},\ldots,x_{n})$, we can
calculate its penalty. The penalty can be mapped to its probability:

$\mathbb{P}(x_{1},\ldots,x_{n})=\frac{e^{-\beta
P(x_{1},\ldots,x_{n})}}{\sum_{x_{1}^{\prime},\ldots,x_{n}^{\prime}}e^{-\beta
P(x_{1}^{\prime},\ldots,x_{n}^{\prime})}},$

where $\beta$ is a parameter that should be determined. We can output the
result configuration with the largest probability, or take expectation for all
configurations to output the probability distribution for each experiment. In
some situations, known results are not deterministic, but stochastic. For
example, assume each experiment has two possible results: $N$ (normal) and $A$
(abnormal). Assume one experiment has result $\mathbb{P}(N)=1/3$, and the
other experiment has result $\mathbb{P}(N)=2/3$. We can assume those two
experiments are independent, and sample deterministic results. Then
$\mathbb{P}(N,N)=2/9$, $\mathbb{P}(N,A)=1/9$, $\mathbb{P}(A,N)=4/9$,
$\mathbb{P}(A,A)=2/9$. For each possible configuration, apply the inference
method, and then take expectations for all configurations.

The above inference method can be applied in the opposite direction. Assume we
want to know the results of some transplantation experiments. We only need to
perform some experiments, and use the inference method to determine other
results. Now the question is to design experiments, so that all the results
can be determined, while the cost is minimized. After transforming the
similarities between experiments into the similarity chart, the experiment
design problem becomes a coloring problem. See Table 2 for the optimal
coloring configurations in two-dimensional lattice, where each not colored
unit (not conducted experiment) is neighboring to $4/2/1$ colored units
(conducted experiment). Since the actual transplantation experiment should be
described by at least four factors (donor type, donor stage, host type, host
stage), the coloring problem should be generalized to any $n$-dimensional
lattice.

|  |  |  |  | Donor |  |
---|---|---|---|---|---|---|---
|  | $T_{1}$ | $T_{2}$ | $T_{3}$ | $T_{4}$ | $T_{5}$ | $T_{6}$ | $T_{7}$
| $T_{1}$ |  |  |  |  |  |  |
| $T_{2}$ |  |  |  |  |  |  |
| $T_{3}$ |  |  |  |  |  |  |
Host | $T_{4}$ |  |  |  |  |  |  |
| $T_{5}$ |  |  |  |  |  |  |
| $T_{6}$ |  |  |  |  |  |  |
| $T_{7}$ |  |  |  |  |  |  |
|  | T1 | T2 | T3 | T4 | T5 | T6 | T7
| T1 |  |  |  |  |  |  |
| T2 |  |  |  |  |  |  |
| T3 |  |  |  |  |  |  |
Host | T4 |  |  |  |  |  |  |
| T5 |  |  |  |  |  |  |
| T6 |  |  |  |  |  |  |
| T7 |  |  |  |  |  |  |
|  | T1 | T2 | T3 | T4 | T5 | T6 | T7
| T1 |  |  |  |  |  |  |
| T2 |  |  |  |  |  |  |
| T3 |  |  |  |  |  |  |
Host | T4 |  |  |  |  |  |  |
| T5 |  |  |  |  |  |  |
| T6 |  |  |  |  |  |  |
| T7 |  |  |  |  |  |  |
Table 2: Experimental designs in different cases, corresponding to two-
dimensional experiment similarity chart $\mathbb{Z}^{2}$. Each cell is an
experiment, and neighboring cells are similar experiments. Black cells are
conducted experiments, and white cells are non-conducted experiments. Black
cells are not neighboring. Each non-boundary white cell is neighboring to $k$
black cells, where $k=4$ (upper), $k=2$ (middle), $k=1$ (lower) [64]. This
table corresponds to the contents of Subsection 4.1.

In this general setting, we need to color each node of the $n$-dimensional
square lattice $\mathbb{Z}^{n}$ black or white, so that each white node is
neighboring to at least $k$ black node, and the number of black nodes is
minimized. Here each node has a coordinate $(x_{1},x_{2},\ldots,x_{n})$, and
two nodes are neighboring if one component of their coordinates differs by
$1$, and other components are equal. Each white node is neighboring to at
least $k$ black nodes, and each black node is neighboring to at most $2n$
white nodes. Therefore, the theoretical upper bound of white-black ratio is
$2n:k$, and the minimal proportion of black nodes is $k/(2n+k)$. The following
results solve the coloring problem for certain $k$ [57]:

######  Theorem 10.

If $k$ is a divisor of $n$, color a node black if and only if its coordinate
satisfies $a_{1}x_{1}+a_{2}x_{2}+\cdots+a_{n}x_{n}\equiv 0\text{ }(\text{mod
}(2n/k)+1)$, where the coefficients $a_{1},a_{2},\ldots,a_{n}$ are $k$ groups
of $1,2,\ldots,n/k$: $1,2,\ldots,n/k,1,2,\ldots,n/k,\ldots,1,2,\ldots,n/k$.
Then black nodes are not neighboring, each white node is neighboring to $k$
black nodes, and the proportion of black nodes is $k/(2n+k)$.

######  Theorem 11.

If $k$ is a divisor of $2n$, but not a divisor of $n$, color a node black if
and only if its coordinate satisfies
$a_{1}x_{1}+a_{2}x_{2}+\cdots+a_{n}x_{n}\equiv 0\text{ }(\text{mod
}(2n/k)+1)$, where the coefficients $a_{1},a_{2},\ldots,a_{n}$ are $k/2$
groups of $1,2,\ldots,2n/k$:
$1,2,\ldots,2n/k,1,2,\ldots,2n/k,\ldots,1,2,\ldots,2n/k$. Then black nodes are
not neighboring, each white node is neighboring to $k$ black nodes, and the
proportion of black nodes is $k/(2n+k)$.

This is another example that cancer research can lead to mathematical
progresses.

###  4.2 Other topics in developmental biology

When transplanting cancer cells into an embryo, the position of
transplantation determines the fates of these cancer cells [30]. This leads to
a central problem in developmental biology: why cells at different positions
of an embryo know where they are, and differentiate accordingly into normal
tissues? The notion of positional information is invented to explain this
problem [69]. Nevertheless, the usages of this term cause confusions, since it
might mean different objects in different studies. In a recent work, the
definition of positional information is clarified with some criteria: it is
position-related; it is extracellular; it is the direct cause of cell fate
[58]. This method of clarifying a notion with criteria can be applied to other
fields [53].

To study the development of an embryo, one convenient way is to represent the
development process as a rooted tree, where each node is a cell along with its
cell event (label), while an edge links the parent cell and the child cell.
Besides, we cannot distinguish two children cells after division. Therefore,
each development process is a rooted unordered tree with possibly repeated
labels. To compare such trees, there are different metrics, where some can be
calculated within polynomial time [54]. One straightforward idea is to
calculate the minimal number of operations to transform one tree into another.
Another idea is to transform both trees into some regular forms and compare
them.

##  5 Discussion

In this review, we go through three facets of mathematical cancer biology:
population dynamics, gene regulation, and developmental biology. These fields
have many meaningful mathematical problems, and solving them can help with the
understanding of cancer. Besides, skills in mathematical modeling and data
analysis are also essential. There are other fields related to cancer biology
that need applied mathematicians. Although most mathematical results in this
review are related to probability and discrete mathematics, other mathematical
tools might be useful in such fields, such as differential equations and
dynamical system, and even geometry. We hope that this review can inspire more
applied mathematicians to work on mathematical problems related to cancer
biology.

## References

  * [1] Altrock, P. M., Liu, L. L., and Michor, F. The mathematics of cancer: integrating quantitative models.  Nature Reviews Cancer 15, 12 (2015), 730–745.
  * [2] Anderson, A. R., Tomlin, C. J., Couch, J., and Gallahan, D. Mathematics of the integrative cancer biology program, 2013.
  * [3] Angelini, E., Wang, Y., Zhou, J. X., Qian, H., and Huang, S. A model for the intrinsic limit of cancer therapy: Duality of treatment-induced cell death and treatment-induced stemness.  PLOS Computational Biology 18, 7 (2022), e1010319.
  * [4] Athreya, K. B., Ney, P. E., and Ney, P. Branching processes.  Courier Corporation, 2004.
  * [5] Ballestar, E., and Esteller, M. Epigenetic gene regulation in cancer.  Advances in genetics 61 (2008), 247–267.
  * [6] Bansal, M., Gatta, G. D., and Di Bernardo, D. Inference of gene regulatory networks and compound mode of action from time course gene expression profiles.  Bioinformatics 22, 7 (2006), 815–822.
  * [7] Bartoszynski, R., Brown, B. W., McBride, C. M., and Thompson, J. R. Some nonparametric techniques for estimating the intensity function of a cancer related nonstationary Poisson process.  The Annals of Statistics (1981), 1050–1060.
  * [8] Biava, P. M., and Bonsignorio, D. Cancer and cell differentiation: a model to explain malignancy.  Journal of Tumor Marker Oncology 17, 2 (2002), 47–54.
  * [9] Bocci, F., Zhou, P., and Nie, Q. spliceJAC: transition genes and state-specific gene regulation from single-cell transcriptome data.  Molecular Systems Biology 18, 11 (2022), e11176.
  * [10] Cao, Z., and Grima, R. Linear mapping approximation of gene regulatory networks with stochastic dynamics.  Nature communications 9, 1 (2018), 1–15.
  * [11] Chaffer, C. L., Marjanovic, N. D., Lee, T., Bell, G., Kleer, C. G., Reinhardt, F., D’Alessio, A. C., Young, R. A., and Weinberg, R. A. Poised chromatin at the zeb1 promoter enables breast cancer cell plasticity and enhances tumorigenicity.  cell 154, 1 (2013), 61–74.
  * [12] Chen, X., and Li, D. ERVcaller: identifying polymorphic endogenous retrovirus and other transposable element insertions using whole-genome sequencing data.  Bioinformatics 35, 20 (2019), 3913–3922.
  * [13] Chen, X., Wang, Y., Feng, T., Yi, M., Zhang, X., and Zhou, D. The overshoot and phenotypic equilibrium in characterizing cancer dynamics of reversible phenotypic plasticity.  Journal of Theoretical Biology 390 (2016), 40–49.
  * [14] Dalerba, P., Cho, R. W., Clarke, M. F., et al. Cancer stem cells: models and concepts.  Annual review of medicine 58, 1 (2007), 267–284.
  * [15] Dessalles, R., Pan, Y., Xia, M., Maestrini, D., D’Orsogna, M. R., and Chou, T. How naive T-cell clone counts are shaped by heterogeneous thymic output and homeostatic proliferation.  Frontiers in immunology 12 (2021).
  * [16] Dewanji, A., Luebeck, E., and Moolgavkar, S. H. A generalized Luria–Delbrück model.  Mathematical biosciences 197, 2 (2005), 140–152.
  * [17] Fessler, E., Borovski, T., and Medema, J. P. Endothelial cells induce cancer stem cell features in differentiated glioblastoma cells via bfgf.  Molecular cancer 14, 1 (2015), 1–12.
  * [18] Gardner, E. J., Lam, V. K., Harris, D. N., Chuang, N. T., Scott, E. C., Pittard, W. S., Mills, R. E., Devine, S. E., Consortium, . G. P., et al. The mobile element locator tool (MELT): population-scale mobile element discovery and biology.  Genome research 27, 11 (2017), 1916–1929.
  * [19] Gatenby, R. A., and Maini, P. K. Mathematical oncology: cancer summed up.  Nature 421, 6921 (2003), 321–321.
  * [20] Glymour, M., Pearl, J., and Jewell, N. P. Causal inference in statistics: A primer.  John Wiley & Sons, 2016.
  * [21] Gupta, P. B., Fillmore, C. M., Jiang, G., Shapira, S. D., Tao, K., Kuperwasser, C., and Lander, E. S. Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells.  Cell 146, 4 (2011), 633–644.
  * [22] Hafner, A., Bulyk, M. L., Jambhekar, A., and Lahav, G. The multiple mechanisms that regulate p53 activity and cell fate.  Nature reviews Molecular cell biology 20, 4 (2019), 199–210.
  * [23] Hurley, D., Araki, H., Tamada, Y., Dunmore, B., Sanders, D., Humphreys, S., Affara, M., Imoto, S., Yasuda, K., Tomiyasu, Y., et al. Gene network inference and visualization tools for biologists: application to new human transcriptome datasets.  Nucleic acids research 40, 6 (2012), 2377–2398.
  * [24] Huynh-Thu, V. A., Irrthum, A., Wehenkel, L., and Geurts, P. Inferring regulatory networks from expression data using tree-based methods.  PloS one 5, 9 (2010), e12776.
  * [25] Jacobs, M., Kim, I., and Tong, J. Tumor growth with nutrients: Regularity and stability.  arXiv preprint arXiv:2204.07572 (2022).
  * [26] Janson, S. Functional limit theorems for multitype branching processes and generalized Pólya urns.  Stochastic Processes and their Applications 110, 2 (2004), 177–245.
  * [27] Jia, C., Qian, H., Chen, M., and Zhang, M. Q. Relaxation rates of gene expression kinetics reveal the feedback signs of autoregulatory gene networks.  The journal of Chemical physics 148, 9 (2018), 095102.
  * [28] Jiang, D.-Q., Wang, Y., and Zhou, D. Phenotypic equilibrium as probabilistic convergence in multi-phenotype cell population dynamics.  PLOS ONE 12, 2 (2017), e0170916.
  * [29] Jordan, C. T., Guzman, M. L., and Noble, M. Cancer stem cells.  New England Journal of Medicine 355, 12 (2006), 1253–1261.
  * [30] Kagawa, H., Javali, A., Khoei, H. H., Sommer, T. M., Sestini, G., Novatchkova, M., Scholte op Reimer, Y., Castel, G., Bruneau, A., Maenhoudt, N., et al. Human blastoids model blastocyst development and implantation.  Nature 601, 7894 (2022), 600–605.
  * [31] Kang, Y., Gu, C., Yuan, L., Wang, Y., Zhu, Y., Li, X., Luo, Q., Xiao, J., Jiang, D., Qian, M., et al. Flexibility and symmetry of prokaryotic genome rearrangement reveal lineage-associated core-gene-defined genome organizational frameworks.  mBio 5 (2014), e01867.
  * [32] Karlin, S. A first course in stochastic processes.  Academic press, 2014.
  * [33] Kent, D. G., and Green, A. R. Order matters: the order of somatic mutations influences cancer evolution.  Cold Spring Harbor perspectives in medicine 7, 4 (2017), a027060.
  * [34] Levine, A. J., Jenkins, N. A., and Copeland, N. G. The roles of initiating truncal mutations in human cancers: the order of mutations and tumor cell type matters.  Cancer Cell 35, 1 (2019), 10–15.
  * [35] Li, H., Ning, S., Ghandi, M., Kryukov, G. V., Gopal, S., Deik, A., Souza, A., Pierce, K., Keskula, P., Hernandez, D., et al. The landscape of cancer cell line metabolism.  Nature medicine 25, 5 (2019), 850–860.
  * [36] Li, W., Wang, Z., Syed, S., Lyu, C., Lincoln, S., O’Neil, J., Nguyen, A. D., Feng, I., and Young, M. W. Chronic social isolation signals starvation and reduces sleep in drosophila.  Nature 597, 7875 (2021), 239–244.
  * [37] Longo, G. Information and causality: Mathematical reflections on cancer biology.  Organisms. Journal of Biological Sciences 2, 1 (2018), 83–104.
  * [38] Medstrand, P., Van de Lagemaat, L., Dunn, C. A., Landry, J.-R., Svenback, D., and Mager, D. L. Impact of transposable elements on the evolution of mammalian gene regulation.  Cytogenetic and genome research 110, 1-4 (2005), 342–352.
  * [39] Meyer, M. J., Fleming, J. M., Ali, M. A., Pesesky, M. W., Ginsburg, E., and Vonderhaar, B. K. Dynamic regulation of cd24 and the invasive, cd44poscd24negphenotype in breast cancer cell lines.  Breast Cancer Research 11, 6 (2009), 1–14.
  * [40] Niu, Y., Wang, Y., and Zhou, D. The phenotypic equilibrium of cancer cells: From average-level stability to path-wise convergence.  Journal of Theoretical Biology 386 (2015), 7–17.
  * [41] Norris, J. R. Markov chains.  Cambridge university press, 1998.
  * [42] Ortmann, C. A., Kent, D. G., Nangalia, J., Silber, Y., Wedge, D. C., Grinfeld, J., Baxter, E. J., Massie, C. E., Papaemmanuil, E., Menon, S., et al. Effect of mutation order on myeloproliferative neoplasms.  New England Journal of Medicine 372, 7 (2015), 601–612.
  * [43] Pollicott, M., Wang, H., and Weiss, H. Extracting the time-dependent transmission rate from infection data via solution of an inverse ode problem.  Journal of biological dynamics 6, 2 (2012), 509–523.
  * [44] Postovit, L.-M., Margaryan, N. V., Seftor, E. A., Kirschmann, D. A., Lipavsky, A., Wheaton, W. W., Abbott, D. E., Seftor, R. E., and Hendrix, M. J. Human embryonic stem cell microenvironment suppresses the tumorigenic phenotype of aggressive cancer cells.  Proceedings of the National Academy of Sciences 105, 11 (2008), 4329–4334.
  * [45] Qian, H., Cheng, Y.-C., and Yang, Y.-J. Kinematic basis of emergent energetics of complex dynamics.  EPL (Europhysics Letters) 131, 5 (2020), 50002.
  * [46] Quintana, E., Shackleton, M., Foster, H. R., Fullen, D. R., Sabel, M. S., Johnson, T. M., and Morrison, S. J. Phenotypic heterogeneity among tumorigenic melanoma cells from patients that is reversible and not hierarchically organized.  Cancer cell 18, 5 (2010), 510–523.
  * [47] Reya, T., Morrison, S. J., Clarke, M. F., and Weissman, I. L. Stem cells, cancer, and cancer stem cells.  nature 414, 6859 (2001), 105–111.
  * [48] Sanchez-Castillo, M., Blanco, D., Tienda-Luna, I. M., Carrion, M., and Huang, Y. A bayesian framework for the inference of gene regulatory networks from time and pseudo-time series data.  Bioinformatics 34, 6 (2018), 964–970.
  * [49] Statello, L., Guo, C.-J., Chen, L.-L., and Huarte, M. Gene regulation by long non-coding RNAs and its biological functions.  Nature reviews Molecular cell biology 22, 2 (2021), 96–118.
  * [50] Tseng, H.-H., Luo, Y., Cui, S., Chien, J.-T., Ten Haken, R. K., and Naqa, I. E. Deep reinforcement learning for automated radiation adaptation in lung cancer.  Medical physics 44, 12 (2017), 6690–6705.
  * [51] Veerman, F., Popović, N., and Marr, C. Parameter inference with analytical propagators for stochastic models of autoregulated gene expression.  International Journal of Nonlinear Sciences and Numerical Simulation 23, 3-4 (2022), 565–577.
  * [52] Wang, Y. Some Problems in Stochastic Dynamics and Statistical Analysis of Single-Cell Biology of Cancer.  Ph.D. thesis, University of Washington, 2018.
  * [53] Wang, Y. Impossibility results about inheritance and order of death.  PLOS ONE 17, 11 (2022), e0277430.
  * [54] Wang, Y. Two metrics on rooted unordered trees with labels.  Algorithms for Molecular Biology 17, 1 (2022), 1–17.
  * [55] Wang, Y. Longest common subsequence algorithms and applications in determining transposable genes.  arXiv preprint arXiv:2301.03827 (2023).
  * [56] Wang, Y., Dessalles, R., and Chou, T. Modelling the impact of birth control policies on China’s population and age: effects of delayed births and minimum birth age constraints.  Royal Society Open Science 9, 6 (2022), 211619.
  * [57] Wang, Y., and He, S. Inference on autoregulation in gene expression.  arXiv preprint arXiv:2201.03164 (2022).
  * [58] Wang, Y., Kropp, J., and Morozova, N. Biological notion of positional information/value in morphogenesis theory.  International Journal of Developmental Biology 64, 10-11-12 (2020), 453–463.
  * [59] Wang, Y., Minarsky, A., Penner, R., Soulé, C., and Morozova, N. Model of morphogenesis.  Journal of Computational Biology 27, 9 (2020), 1373–1383.
  * [60] Wang, Y., Mistry, B. A., and Chou, T. Discrete stochastic models of selex: Aptamer capture probabilities and protocol optimization.  The Journal of Chemical Physics 156, 24 (2022), 244103.
  * [61] Wang, Y., and Qian, H. Mathematical representation of Clausius’ and Kelvin’s statements of the second law and irreversibility.  Journal of Statistical Physics 179, 3 (2020), 808–837.
  * [62] Wang, Y., and Wang, L. Causal inference in degenerate systems: An impossibility result.  In International Conference on Artificial Intelligence and Statistics (2020), PMLR, pp. 3383–3392.
  * [63] Wang, Y., and Wang, Z. Inference on the structure of gene regulatory networks.  Journal of Theoretical Biology 539 (2022), 111055.
  * [64] Wang, Y., Zhang, B., Kropp, J., and Morozova, N. Inference on tissue transplantation experiments.  Journal of Theoretical Biology 520 (2021), 110645.
  * [65] Wang, Y., and Zheng, Z. Measuring policy performance in online pricing with offline data.  Available at SSRN 3729003 (2021).
  * [66] Wang, Y., Zheng, Z., and Shen, Z.-J. M. Online pricing with polluted offline data.  Available at SSRN 4320324 (2023).
  * [67] Wang, Y., Zhou, J. X., Pedrini, E., Rubin, I., Khalil, M., Qian, H., and Huang, S. Multiple phenotypes in HL60 leukemia cell population.  arXiv preprint arXiv:2301.03782 (2023).
  * [68] Wodarz, D., and Komarova, N. Computational biology of cancer: lecture notes and mathematical modeling.  World Scientific, 2005.
  * [69] Wolpert, L. Positional information and pattern formation.  Current topics in developmental biology 117 (2016), 597–608.
  * [70] Xia, M., Li, X., and Chou, T. An age-structured Lotka-Volterra model and the emergence of overcompensation.  Bulletin of the American Physical Society (2023).
  * [71] Xing, B., and Van Der Laan, M. J. A causal inference approach for constructing transcriptional regulatory networks.  Bioinformatics 21, 21 (2005), 4007–4013.
  * [72] Yang, G., Quan, Y., Wang, W., Fu, Q., Wu, J., Mei, T., Li, J., Tang, Y., Luo, C., Ouyang, Q., et al. Dynamic equilibrium between cancer stem cells and non-stem cancer cells in human sw620 and mcf-7 cancer cell populations.  British journal of cancer 106, 9 (2012), 1512–1519.
  * [73] Yang, Y.-J., and Cheng, Y.-C. Potentials of continuous Markov processes and random perturbations.  Journal of Physics A: Mathematical and Theoretical 54, 19 (2021), 195001.
  * [74] Ye, F. X.-F., Wang, Y., and Qian, H. Stochastic dynamics: Markov chains and random transformations.  Discrete & Continuous Dynamical Systems-B 21, 7 (2016), 2337.
  * [75] Yu, T., Huang, X., Dou, S., Tang, X., Luo, S., Theurkauf, W. E., Lu, J., and Weng, Z. A benchmark and an algorithm for detecting germline transposon insertions and measuring de novo transposon insertion frequencies.  Nucleic Acids Research 49, 8 (2021), e44–e44.
  * [76] Zhao, Y., Kosorok, M. R., and Zeng, D. Reinforcement learning design for cancer clinical trials.  Statistics in medicine 28, 26 (2009), 3294–3315.
  * [77] Zhao, Y., Zeng, D., Socinski, M. A., and Kosorok, M. R. Reinforcement learning strategies for clinical trials in nonsmall cell lung cancer.  Biometrics 67, 4 (2011), 1422–1433.
  * [78] Zhou, D., Wang, Y., and Wu, B. A multi-phenotypic cancer model with cell plasticity.  Journal of Theoretical Biology 357 (2014), 35–45.
  * [79] Zhou, P., Wang, S., Li, T., and Nie, Q. Dissecting transition cells from single-cell transcriptome data through multiscale stochastic dynamics.  Nature communications 12, 1 (2021), 1–15.


# On the Nonlinear Theory of Viscoelasticity of Differential Type

Edvige Pucci, Giuseppe Saccomandi,
Dipartimento di Ingegneria Industriale,
Università degli Studi di Perugia, 06125 Perugia, Italy.

###### Abstract

We consider nonlinear viscoelastic materials of differential type and for some
special models we derive exact solutions of initial boundary value problems.
These exact solutions are used to investigate the reasons of non-existence of
global solutions for such equations.

Keywords: Viscoelasticity of differential type; nonlinear viscosity functions,
exact solutions, existence and uniqueness of boundary value problems.

##  1 Introduction and basic equations

In empirical rheology the standard sketch of the one dimensional Kelvin Voigt
model of linear viscoelasticity is of a purely viscous damper and a purely
elastic spring connected in parallel. In continuum mechanics, this model may
be recast in the class of the simple materials denoted as viscoelastic solids
of differential type. For such materials, in the isotropic case, the Cauchy
stress tensor, $\boldsymbol{T}$, is defined by the constitutive equation

$\boldsymbol{T}=\boldsymbol{T}(\boldsymbol{B};\boldsymbol{D}),$ (1.1)

where $\boldsymbol{B}=\boldsymbol{F}\boldsymbol{F}^{T}$ and
$\boldsymbol{D}=(\dot{\boldsymbol{F}}\boldsymbol{F}^{-1}+\boldsymbol{F}^{-T}\dot{\boldsymbol{F}}^{T})/2$
are the left Cauchy-Green strain tensor and the stretching tensor respectively
(Truesdell and Noll, 1965). (Here $\boldsymbol{F}$ is the gradient of
deformation). Therefore $\boldsymbol{B}$ stands for the elastic spring and
$\boldsymbol{D}$ stands for the viscous damper.

Nowadays there is an increasing interest in the study of nonlinear models of
viscoelasticity and this because in the experimental studies of rubber-like
materials and soft tissues the viscous and nonlinear effects are prominent and
clearly measurable (Lakes, 2009). For these reasons the possibility to
elaborate feasible and reliable mathematical models of nonlinear
viscoelasticity is a major problem of contemporary continuum mechanics
(Wineman, 2009). In this framework the class of constitutive equations (1.1)
seems to have a central role in view of its direct and clear connection with
the linear Kelvin-Voigt model.

The aim of the present note is to show that many of the nonlinear models
contained in the constitutive equation (1.1) are not _acceptable_ because they
are clearly unphysical.

To make the algebra simple as possible we shall restrict our attention to the
following family of shearing motions

$x=X+u(Z,t),\qquad y=Y,\qquad z=Z.$ (1.2)

Moreover, it is convenient to consider only incompressible materials. In the
case of compressible materials some complications arise, because the family of
motions (1.2) is in general not admissible (Destrade and Saccomandi, 2005) and
because the constitutive issues are more complex.

In this framework, it is straightforward to check that the balance equations
(in the absence of body forces) reduce to the scalar equations

$\rho u_{tt}=\partial T_{13}/\partial Z+p_{x},\quad p_{y}=0,\quad p_{z}=0,$
(1.3)

where $\rho$ is the (constant) density of the material and $p$ is the Lagrange
multiplier associated with the constraint of incompressibility
($det\boldsymbol{F}=1$). For an isotropic material the classical
representation formula (Truesdell and Noll, 1965) gives

$T_{13}=Q(K^{2})+\nu(K^{2},K_{t}^{2},KK_{t})K_{t}$ (1.4)

where $K(Z,t)=u_{Z}$ is the amount of shear, $Q=Q(K^{2})$ is the generalized
shear modulus and $\nu=\nu(K^{2},K_{t}^{2},KK_{t})$ is the viscosity
_function_. The unknowns of our problem are the pressure field $p=p(x,y,z,t)$
and the shearing motion $u=u(Z,t)$ (or equivalently the amount of shear $K$).

From the three scalar equations (1.3) we first deduce that

$p(x,t)=f(t)x+g(t),$

where $g(t)$ is an unessential function. Then by derivation, with respect of
$Z$, of the first equation in (1.3) we obtain

$\rho K_{tt}=\frac{\partial}{\partial
Z^{2}}\left[Q(K^{2})K+\nu(K^{2},K_{t}^{2},KK_{t})K_{t}\right].$ (1.5)

Once that (1.5) has been solved for $K(Z,t)$ and we integrate to obtain
$u(Z,t)$ from (1.3) it is possible to recover the pressure field, i.e. the
function $f(t)$.

Assuming that there is a characteristic time, $\omega^{-1}$, such that
$\tau=\omega t$, a characteristic length, $L$ ,such that $\zeta=Z/L$ and
assuming that the limit

$\mu_{0}=lim_{K\rightarrow 0}Q$

is a positive number, we define the quantities

$\epsilon=\frac{\omega^{2}L^{2}\rho}{\mu_{0}},\quad\hat{Q}=\frac{Q}{\mu_{0}},\quad\hat{\nu}=\frac{\omega\nu}{\mu_{0}},$

such that the dimensionless version of (1.5) reads

$\epsilon
K_{\tau\tau}=\frac{\partial}{\partial\zeta^{2}}\left[\hat{Q}(K^{2})K+\hat{\nu}(K^{2},\omega^{2}K_{\tau}^{2},\omega
KK_{\tau})K_{\tau}\right].$ (1.6)

When $\hat{Q}=1$, $\nu=\nu_{0}$ choosing $\omega=\mu_{0}/\nu_{0}$, (i.e.
$\hat{\nu}=1$) from (1.6) we recover the classical _linear_ Kelvin-Voigt
material

$\epsilon K_{\tau\tau}=K_{\zeta\zeta}+K_{\tau\zeta\zeta}..$ (1.7)

If in (1.7) we let $\epsilon\rightarrow 0$ we obtain the _quasistatic_
approximation. This approximation is fundamental to model the usual creep and
recovery experimental tests (Lakes, 2009). In Pucci and Saccomandi (2010) it
has been shown that when:

  * •

the empirical inequalities are in force111In this case the generalized shear
modulus, $Q$, is positive for any value of the amount of shear and the strain-
energy density function associated with the elastic part of equation (1.6)
must be convex. (Beatty, 1987);

  * •

the infinitesimal shear modulus $\mu_{0}$ is positive and finite;

  * •

the viscosity function is such that

$\lim_{K\rightarrow 0,K_{t}\rightarrow 0}\nu=\nu_{0},$

where $\nu_{0}$ is a positive constant;

then the constitutive equations (1.1) are able to describe, at least
qualitatively, the mechanical behavior observed in standard quasi-static
experimental tests.

The quasi-static approximation is a _singular perturbation_ of the full
dynamical equation (1.6) and therefore is an approximation valid for large
times. For this reason the global existence of solutions of (1.6) is a main
issue. The existence, uniqueness and regularity of solutions for initial-
boundary value problems (IBVPs) associated to (1.5) (or its dimensionless
version (1.6)) have been studied by several authors. A state of art of the
various results both for classical and weak solutions is provided in (Tvedt,
2008).

The main problem with (1.5) is that the existence of smooth solutions global
in time is not guaranteed when the viscosity function is nonlinear. Indeed
Tvedt points out that:

> To my knowledge, all existence results for weak solutions of any kind to
> (1.1) in multiple space dimensions, including the Young-measure valued
> constructed in [6], require at least
>
> $\left|S(p,q)\right|\leq C(1+\left|p\right|+\left|q\right|)$

> for some constant $C$. _This condition is difficult to justify on physical
> grounds_.

(Using the present notation $S=T_{13}$, $p=K$, $q=K_{t}$.).

It seems that the consequences of this negative result quoted by Tvedt are not
full appreciated by workers in nonlinear viscoelasticity. For example, if we
have not a global existence result it is clear that we can consider a quasi-
static approximation. Our show note is centered on the impact of this result
in very simple problems of viscoelasticity.

The plan of the paper is the following. In the next Section we derive some
simple exact solutions for simple IBVPs with very smooth initial and boundary
conditions of special cases of (1.6) . For all these solutions there is a
finite time $t^{*}$ associated with the blow-up or the extinction of the
solution. In Section 3, we suggest that the _catastrophic_ problem associated
with $t^{*}$ is located in the initial layer of the solution (i.e. for short
times). By considering the case where the viscosity function depends only on
the amount of shear we prove that the dynamics in the initial layer is
governed by an equation equivalent to a nonlinear diffusion equation. The
blow-up characteristics of such an equation have been studied into details in
(Galaktionov and Vazquez, 1995) and the conditions derived there show the non-
existence of global solutions. The last Section is devoted to concluding
remarks.

##  2 Some exact solutions

We start by considering the problem associated with some viscous function that
does not have an upper bound. To this end we consider the following three
nonlinear polynomial viscous functions in turn

$\hat{\nu}_{e}(K^{2})=1+\hat{\nu}_{e1}K^{2},\quad\hat{\nu}_{v}(K_{\tau}^{2})=1+\hat{\nu}_{v1}K_{\tau}^{2},\quad\hat{\nu}_{m}(KK_{\tau})=1+\hat{\nu}_{m1}KK_{\tau}.$

For all these viscosity functions we have that $\lim_{K\rightarrow
0,K_{\tau}\rightarrow 0}\hat{\nu}=1>0$ but none of these functions are bounded
from above.

We start considering, for $\zeta\in[0,1]$, the equation

$\epsilon
K_{\tau\tau}=K_{\zeta\zeta}+K_{\tau\zeta\zeta}+\hat{\nu}_{e1}\left(K^{2}K_{\tau}\right)_{\zeta\zeta},$
(2.1)

i.e. an equation corresponding to a material whose elastic part is neo-Hookean
and with a super linear viscosity function in the amount of shear defined as
$\hat{\nu}_{e}$.

This equation has an interesting and simple exact solution that may be
computed as

$K(\zeta,\tau)=\alpha(\tau)\zeta+\beta(\tau),$ (2.2)

where

$\epsilon\frac{d\alpha}{d\tau}=2\hat{\nu}_{e1}\alpha^{3}+c_{1},\quad\epsilon\frac{d\beta}{d\tau}=2\hat{\nu}_{e1}\alpha^{2}\beta+c_{2},$

(here $c_{1}$ and $c_{2}$ are two constants of integration). If we consider
$c_{1}=c_{2}=0$ and

$\beta(\tau)\equiv-\frac{1}{2}\alpha(\tau)\equiv\frac{1}{2}\sqrt{\frac{\epsilon}{\epsilon-\hat{\nu}_{e1}\tau}}$

we obtain

$K(\zeta,\tau)=\frac{1}{2}\sqrt{\frac{\epsilon}{\epsilon-\hat{\nu}_{e1}\tau}}(\zeta-\frac{1}{2})\rightarrow\hat{u}(\zeta,\tau)=\frac{1}{4}\sqrt{\frac{\epsilon}{\epsilon-\hat{\nu}_{e1}\tau}}(\zeta^{2}-\zeta).$
(2.3)

(Here $\hat{u}=u/L$).

The shearing motion in (2.3) satisfies the following set of smooth initial and
boundary conditions

$\hat{u}(0,\tau)=\hat{u}(1,\tau)=0,\,\hat{u}(\zeta,0)=(\zeta^{2}-\zeta)/4,\,\hat{u}_{\tau}(\zeta,0)=\hat{\nu}_{1}(\zeta^{2}-\zeta)/(8\epsilon).$

In this way we have provided a simple example of non-existence (because of
finite time blow-up) of the solution for a very regular IVBPs for equation
(2.1) in $[0,1]$.

If now we consider the equation

$\epsilon
K_{\tau\tau}=K_{\zeta\zeta}+K_{\tau\zeta\zeta}+\hat{\nu}_{v1}\left(K_{\tau}^{3}\right)_{\zeta\zeta},$
(2.4)

a similar solution may be obtained using the ansatz in (2.2) for the amount of
shear. The exact solution is now

$\hat{u}(\zeta,\tau)=\frac{1}{2}\sqrt{\frac{12\hat{\nu}_{v1}-\epsilon\tau}{3\hat{\nu}_{v1}}}(\zeta^{2}-\zeta).$
(2.5)

This solution may be used to solve the boundary-initial value problem

$\hat{u}(0,\tau)=\hat{u}(1,\tau)=0,\,\hat{u}(\zeta,0)=\zeta^{2}-\zeta,\,\hat{u}_{\tau}(\zeta,0)=-\frac{\epsilon}{24\hat{\nu}_{{v1}}}(\zeta^{2}-\zeta).$

Here we record that the reason why the existence of the solution fails in
finite time is the extinction phenomena.

The last equation we consider is

$\epsilon
K_{\tau\tau}=K_{\zeta\zeta}+K_{\tau\zeta\zeta}+\hat{\nu}_{m1}\left(KK_{\tau}^{2}\right)_{\zeta\zeta}.$
(2.6)

Here the mathematics is more involved. Let us consider that the shearing
motion is in the form

$\hat{u}(\zeta,\tau)=\alpha(\tau)(\zeta^{2}-\zeta).$

Then the time dependent function $\alpha$ is determined by the ordinary
differential equation

$\frac{\epsilon}{\hat{\nu}_{m1}}\alpha^{\prime\prime}=24\alpha\alpha^{\prime
2}.$

The solution of this second order ordinary differential equation is given by

$\text{Erf}\left(\tilde{\alpha}\right)=c_{1}\tau+c_{2},$

where

$\tilde{\alpha}=2\sqrt{\frac{3\hat{\nu}_{m1}}{\epsilon}}\alpha,$

Erf is the standard error function and $c_{1},\,c_{2}$ two integration
constants. This exact solution solves the following IBVP (in the normalized
variables)

$\tilde{u}(0,\tau)=\tilde{u}(1,\tau)=0,\,\tilde{u}(\zeta,0)=\delta(\zeta^{2}-\zeta),$

with $0<\delta<1$ and $\tilde{u}_{\tau}(\zeta,0)=\mu(\zeta^{2}-\zeta)$ where
$\mu$ is an ad hoc constant determined by $\alpha^{\prime}$. From

$\text{Erf}\left(\tilde{\alpha}\right)-\text{Erf}\left(\delta\right)=c_{1}\tau,$
(2.7)

it is clear that this solution blows up in a finite time.

To summarize, the following remarks apply to all the three examples we have
proposed. First of all in the case $\epsilon\rightarrow 0$ all the exact
solutions we propose degenerate and are of no interest. Second, for all three
solutions there exists a critical time $\tau^{*}$ such that for
$\tau>\tau^{*}$ the solution does not exists anymore. For the solution in
(2.3) and (2.7) we record a blow-up phenomenon and for the solution (2.5) we
have an extinction phenomenon.

Therefore we have provided explicit examples of the lack of global existence
for equations with nonlinear viscous function and a neo-Hookean elastic term.
For all these solutions the inertia term (i.e. $\epsilon\neq 0$) is crucial in
delineating why the global existence fails.

##  3 The initial layer

We point out once again that in the quasi-static limit all the models we have
considered in the previous Section are ”good” models. This fact suggests that
because the quasi-static limit is the result of a singular perturbation
procedure, the failure of well-posedness must occur for short times or, using
the terminology of singular perturbations in the _initial layer_.

When we consider a general viscosity function defined as
$\hat{\nu}=\hat{\nu}(K^{2})$ and we introduce the _stretched_ time
$\tilde{\tau}=\tau/\epsilon$, it is possible to obtain easily, for any regular
viscosity function, the equation governing the dynamics in the initial layer.
Restricting the attention to this class of materials it is then possible to
rewrite (1.6) as

$K_{\tilde{\tau}\tilde{\tau}}=\frac{\partial}{\partial\zeta^{2}}\left[\epsilon\hat{Q}(K^{2})K+\hat{\nu}(K^{2})K_{\tilde{\tau}}\right],$
(3.1)

and when $\epsilon\rightarrow 0$ we derive the equation

$K_{\tilde{\tau}\tilde{\tau}}=\frac{\partial^{3}\Pi}{\partial\zeta^{2}\partial\tilde{\tau}}.$
(3.2)

where $\Pi=\int\hat{\nu}(K^{2})dK$.

By a simple integration the above equation is rewritten as

$K_{\tilde{\tau}}=\frac{\partial^{2}\Pi}{\partial\zeta^{2}}+\varphi(\zeta).$
(3.3)

Equation (3.3) is a nonlinear diffusion equation. This equation has been
considered by many authors and the necessary and sufficient conditions for the
complete blow-up and extinction of this equation have been considered in
(Galaktionov and Vazquez, 1995). It is interesting to note that for equation
(3.3) we have blow-up of the solution in finite time if the integral

$\int_{1}^{\infty}\frac{\Pi^{\prime}(s)}{s}ds,$

does not converge. (Here a dash denotes the derivative with respect the amount
of shear). This is exactly the case of the examples in the previous Section.

It is clear that there are several choices of the viscosity functions for
which the above integral may converge, but the usual higher order polynomial
models that are used in the standard fitting procedures do blow up. We remark
that a model with a constant viscosity, in the initial layer, is associated
with a linear heat equation this explain why it is possible to find a unique
global solution for such a model in the class (1.1).

We remark that via a multiple scale methods is possible to deduce the equation
in the initial layer for more general viscosity functions, but the case we
have just examined is sufficient to point out the problems inherent with the
constitutive equation (1.1) when the sub linear growth condition pointed out
by Tvedt (2009) is not satisfied.

##  4 Concluding remarks

If we believe that materials of differential type have a ”good” quasi-static
behavior, the fact that a global solution for the full dynamical equations
does not exist is a major problem. Indeed, if a solution blows up, or is
associated with an extinction phenomena, in finite time it is not possible to
consider the asymptotic limit for $t\rightarrow\infty$ and therefore it is not
possible to consider the quasi-static approximation.

Clearly such situation may be arranged considering a sort of ”cut-off” for the
usual polynomial constitutive equations. For example we may be consider a
constitutive equation such as

$\hat{\nu}(K^{2})=\left\\{\begin{array}[]{ll}1+\nu_{e1}K^{2}&\mbox{if
$K^{2}\leq K^{*}$};\\\ \nu_{e}&\mbox{if $K^{2}>K^{*}$}.\end{array}\right.$

We cannot be satisfied with this approach for two reasons. First of all we
introduce some oddities in the quasi-static regime as it is possible to check
by analyzing the corresponding creep problem that loose regularity. Moreover,
we do not solve the ”unphysical” situation generated by the fact: greater
viscous effects lead to greater problems with the regularity of the solution.

This situation seems to suggest that there is something missing in the
”dynamics” of the nonlinear viscoelastic models we are considering. It seems
that the model (1.1) must be refined taking into account phenomena that at
first sight seemed to be irrelevant but indeed are fundamental. These
phenomena seems to be of the same order of the inertia terms and therefore
disappear when $\epsilon\rightarrow 0$.

A possibility is to consider the presence of objective derivatives of the
stretching tensor $\boldsymbol{D}$ i.e. higher order models of differential
type where the constitutive equations not only depends on the stretching
tensor $\boldsymbol{D}$ but for example also on the subsequent Rivlin-Ericksen
tensors (Truesdell and Noll, 1965). These tensors may introduce some
additional terms that regularize the equations and allows to avoid the
problems presented here. This direction is outlined in the context of non
dissipative materials in (Destrade and Saccomandi, 2006).

## Acknowledgements

The research is supported by PRIN-2009 project _Matematica e meccanica dei
sistemi biologici e dei tessuti molli_.

GS thanks Michel Destrade and Jeremiah Murphy for their remarks on a previous
draft of this paper.

This paper is dedicated with great admiration to Professor Robin Knops for his
80-th birthday.

## References

  * [1] Beatty, M.F. : Topics in finite elasticity: Hyperelasticity of rubber, elastomers, and biological tissues-with examples, Appl. Mech. Rev. 40 (1987) 1699–1734.
  * [2] Destrade, M. and Saccomandi, G.: On finite amplitude elastic waves propagating in compressible solids, Physical Review E, 72 (2005) 016620.
  * [3] Destrade, M. and Saccomandi, G.: Solitary and compact-like shear waves in the bulk of solids, Physical Review E, 73 (2006) 065604.
  * [4] Galaktionov, V.A. and Vazquez, J. L.: Necessary and sufficient conditions for complete blow-up and extinction for one-dimensional quasilinear heat equations. Arch. Rational Mech. Anal. 129 (1995) 225–244.
  * [5] Lakes, R. S.: Viscoelastic Materials, Cambridge University Press, 2009
  * [6] Pucci, E. and Saccomandi, G.: On a Special Class of Nonlinear Viscoelastic Solids. Mathematics and Mechanics of Solids 15 (2010) 803–811.
  * [7] C. Truesdell and W. Noll.: The Non-linear Field Theories of Mechanics. Volume III/3, Encyclopedia of Physics (S. Flugge Ed.). Springer, Berlin, 1965.
  * [8] B. Tvedt.: Quasilinear Equations for Viscoelasticity of Strain-Rate Type. Archive for Rational Mechanics and Analysis 189 (2008) 237–281.
  * [9] Wineman, A.: Nonlinear viscoelastic solids-a review, Mathematics and Mechanics of Solids 14 (2009) 300-366.