Open Quasispecies Models: Stability, Optimization, and Distributed Extension

TL;DR

This paper introduces open quasispecies models that incorporate extinction, analyzing their stability and dynamics both analytically and numerically, and extending them to distributed models for complex biological applications.

Contribution

It proposes a novel open quasispecies framework with stability analysis, inverse fitness optimization, and continuous distributed models for large genotype spaces.

Findings

Quasispecies distribution and error threshold are observed in the models.

Extinction phenomena are incorporated alongside growth dynamics.

Distributed models enable analysis of high-dimensional genotype spaces.

Abstract

We suggest a natural approach that leads to a modification of classical quasispecies models and incorporates the possibility of population extinction in addition to growth. The resulting modified models are called open. Their essential properties, regarding in particular equilibrium behavior, are investigated both analytically and numerically. The hallmarks of the quasispecies dynamics, viz. the heterogeneous quasispecies distribution itself and the error threshold phenomenon, can be observed in our models, along with extinction. In order to demonstrate the flexibility of the introduced framework, we study the inverse problem of fitness allocation under the biologically motivated criterion of steady-state fitness maximization. Having in mind the complexity of numerical investigation of high-dimensional quasispecies problems and the fact that the actual number of genotypes or alleles involved in a studied process can be extremely large, we also build continuous-time distributed open quasispecies models. The obtained results may serve as an initial step to developing mathematical models that involve directed therapy against various pathogens.

Equations

$$\dot{p}(t)=KMp(t)-f[p(t)]p(t),t⩾0,$$

$$\sum _{i=1}^n{k}_{ij}=1,j=\overline{1,n};$$

$$\dot{p}(t)=\left(M+\hat{M}\right)p(t)-f[p(t)]p(t),t⩾0,$$

$${\mu }_{jj}=-\sum _{i\in \{1,2,\dots ,n\}\setminus \{j\}}{\mu }_{ij},$$

$$\hat{M}=\mu Q,\mu >0,Q\in {\mathbb{R}}^{(N+1)\times (N+1)},$$

$$Q=\left(\begin{array}{cccccccc}-N& 1& 0& 0& \dots & \dots & 0& 0\\ N& -N& 2& 0& \dots & \dots & 0& 0\\ 0& N-1& -N& 3& \dots & \dots & 0& 0\\ 0& 0& N-2& -N& \dots & \dots & 0& 0\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0& 0& \dots & \dots & \dots & -N& N-1& 0\\ 0& 0& \dots & \dots & \dots & 2& -N& N\\ 0& 0& \dots & \dots & \dots & 0& 1& -N\end{array}\right).$$

$$\dot{p}(t)=Gp(t)-f[p(t)]p(t),t⩾0,$$

$$\sum _{i=1}^n{g}_{ij}=m_j,j=\overline{1,n},$$

$$g_{ij}⩾0\mathrm{for}i\ne j,i=\overline{1,n},j=\overline{1,n}.$$

$$p(t)=\frac{1}{{\sum }_{i=1}^n{r}_i(t)}r(t)\forall t⩾0,$$

$$r(t)=(r_1(t),r_2(t),\dots ,r_n(t){)}^{\top }=c\exp \left\{\underset{0}{\overset{t}{\int }}f[p(\tau )]\mathrm{d}\tau \right\}p(t)$$

$$\forall t⩾0$$

$$\dot{r}(t)=Gr(t),t⩾0,$$

$$r(0)=cp(0),$$

$$\dot{u}(t)=\phi \left(\sum _{i=1}^n{u}_i(t)\right)Gu(t)-Du(t),t⩾0.$$

$$u(0)=u^0={\left(u_1^0,u_2^0,\dots ,u_n^0\right)}^{\top }\in [0,+\infty {)}^{n\times 1}.$$

$$F_i(u)=\phi \left(\sum _{j=1}^n{u}_j\right)\sum _{j=1}^n{g}_{ij}{u}_j-d_i{u}_i$$

$$\forall u=(u_1,u_2,\dots ,u_n{)}^{\top }\in {\mathbb{R}}^{n\times 1},i=\overline{1,n}.$$

$$\mathrm{D}F(u)={\left\{\frac{\partial F_i(u)}{\partial u_j}\right\}}_{i,j=1}^n,$$

$$\frac{\partial F_i(u)}{\partial u_j}=\{\begin{array}{ll}{\phi }^{\prime }\left(\sum _{\nu =1}^n{u}_{\nu }\right)\sum _{\nu =1}^n{g}_{i\nu }{u}_{\nu }+\phi \left(\sum _{\nu =1}^n{u}_{\nu }\right)g_{ii}-d_i,& i=j,\\ {\phi }^{\prime }\left(\sum _{\nu =1}^n{u}_{\nu }\right)\sum _{\nu =1}^n{g}_{i\nu }{u}_{\nu }+\phi \left(\sum _{\nu =1}^n{u}_{\nu }\right)g_{ij},& i\ne j,\end{array}$$

$$\forall u=(u_1,u_2,\dots ,u_n{)}^{\top }\in {\mathbb{R}}^{n\times 1},i=\overline{1,n},j=\overline{1,n}.$$

$$\phi (s)=e^{-\gamma s}\forall s\in \mathbb{R},\gamma =\mathrm{const}>0.$$

$$s(t)=\sum _{i=1}^n{u}_i(t)\forall t⩾0$$

$$D^{-1}=\mathrm{diag}\left[d_1^{-1},d_2^{-1},\dots ,d_n^{-1}\right].$$

$$\phi (0)=1,l_{\phi }=0,{\phi }^{-1}(\rho )=\frac{1}{\gamma }\ln \frac{1}{\rho }⩾0\forall \rho \in (0,1],$$

$$D^{-1}G{u}^{\ast }=\frac{1}{\phi (s^{\ast })}{u}^{\ast },s^{\ast }=\sum _{i=1}^n{u}_i^{\ast }>0.$$

$$\phi (s^{\ast })=\frac{1}{{\lambda }^{\ast }}\in (l_{\phi },\phi (0)),s^{\ast }=\sum _{i=1}^n{u}_i^{\ast }={\phi }^{-1}\left(\frac{1}{{\lambda }^{\ast }}\right)>0.$$

$${\lambda }^{\ast }=\frac{{\sum }_{i=1}^n{m}_i{u}_i^{\ast }}{{\sum }_{i=1}^n{d}_i{u}_i^{\ast }}.$$

$$\tilde{f}[u]=\{\begin{array}{ll}0,& {\sum }_{i=1}^n{u}_i=0,\\ \frac{f[u]}{{\sum }_{i=1}^n{d}_i{u}_i}=\frac{{\sum }_{i=1}^n{m}_i{u}_i}{{\sum }_{i=1}^n{d}_i{u}_i},& {\sum }_{i=1}^n{u}_i>0,\end{array}$$

$$\forall u=(u_1,u_2,\dots ,u_n{)}^{\top }\in [0,+\infty {)}^{n\times 1}.$$

$${\lambda }^{\ast }=\tilde{f}[u^{\ast }].$$

$${\lambda }^{\ast }>0,\frac{1}{{\lambda }^{\ast }}\in (l_{\phi },\phi (0)).$$

$$s^{\ast }=\sum _{i=1}^n{u}_i^{\ast }={\phi }^{-1}\left(\frac{1}{{\lambda }^{\ast }}\right).$$

Full text

1 Introduction

Mathematical modeling has long been a key theoretical tool connecting various pictorial and verbal models of RNA virus evolution, establishment, and extinction, as well as in vitro and in vivo experiments aimed to understand and potentially fight this rich group of human pathogens [1]. Arguably, the onset of this specific modeling, which was done originally in terms of self-replicating polynucleotide sequences to study the problem of the origin of life, can be traced back to the founding papers by M. Eigen et al. [2, 3, 4, 5, 6], as well as by J. F. Crow and M. Kimura [7, 8, 9]. They proposed so-called quasispecies models, which were later connected directly to RNA virus evolution [10]; see [11, 12] for a more recent account.

For such quasispecies models, two key phenomena were observed. The first one was the highly heterogeneous cloud of mutants of the most fit (or master) sequence at the selection-mutation equilibrium. This cloud was called the quasispecies, hence the name of the model. The second phenomenon was the so-called error threshold. It can be described as the critical mutation rate (or probability, depending on the settings of a particular model) that divides the selective phase of the virus evolution, i. e., the dominance of the master sequence in the population, from the random phase, where the distribution of various genotypes becomes essentially uniform.

A great deal of mathematical investigation was devoted to study exact properties of the quasispecies and error threshold; see, e. g., the review in [13]. In a nutshell, the exact details of the structure of the mutation-selection equilibrium and the precise position (if it exists at all) of the error threshold depend in a subtle way on the implemented fitness landscape. They are in general far from the oversimplified and widely referred formula that the error threshold is typically given by the selective advantage of the master sequence over the sequence length; see, e. g., [14, 15, 16] for specific examples.

From a more practical point of view, the error threshold concept led to the idea of lethal mutagenesis, i. e., the process of virus extinction induced by elevated mutation rates [17], and eventually to first mutagenetic experiments [18]. The undeniable influence of the error threshold phenomena on the idea of lethal mutagenesis notwithstanding, there existed an internal contradiction between these two concepts from the modeling point of view. The fact is that the quasispecies models were formulated as systems of ordinary differential equations (ODEs) describing the distribution of the relative frequencies of the master sequence and its various mutants, but not their population sizes (or densities with respect to a certain external measure). It was therefore simply meaningless to discuss possible extinction within the framework of these models [19, 20]. In [19], it was reasonably noticed that the process of the loss of the master sequence in the population (the error threshold) should in general be distinguished from the actual extinction of a whole virus population, and the original quasispecies models would have to be modified to incorporate the possibility of population extinction. The paper [19] was the first to suggest such a modification. Under a number of simplifying assumptions (e. g., the multiplicativity of the fitness landscape, no back or compensatory mutations, etc.), the following simple mathematical condition for population extinction was obtained: the product of the absolute fitness of the master sequence and the average biological fitness should be less than one; see also [21].

Other examples of quasispecies models including the possibilities of lethal mutations and population extinction were given in [22, 23, 24]. They were reviewed and put in a general context in [21]. Among other things, the analysis (analytical and/or numerical) of these models showed, quite naturally, that taking into account one or another more realistic assumption would lead to corrections of the original extinction criterion obtained in [19]. Moreover, analytical results are usually possible only under some significant simplifications regarding the fitness and mutational landscapes.

The main goal of the present text is to develop a relatively simple but yet general and flexible modeling framework incorporating growth and mortality characteristics, as well as to investigate its basic properties. Due to the possibility of population extinction, our models can be called open, in contrast to the classical quasispecies models (formulated in terms of relative frequencies) which can be called closed. Our approach differs from those presented in the aforementioned works [19, 22, 23, 24, 21], since we do not start with any radical simplifying assumptions, e. g., we do not fix the fitness or mutational landscapes. Our requirement is to stay connected with the classical quasispecies models as close as reasonably possible. Among other things, such a requirement potentially yields the opportunity to support our study with a number of existing mathematical techniques for the classical quasispecies models; see, e. g., [25, 26, 27, 28, 15].

In order to demonstrate the flexibility of our framework in case of ODE-based dynamics, we also study the inverse problem of fitness allocation under the biologically motivated criterion of steady-state fitness maximization.

Furthermore, having in mind the complexity of numerical investigation of high-dimensional quasispecies problems and the fact that the actual number of genotypes or alleles involved in a studied process can be extremely large, we propose continuous-time distributed open quasispecies models with growth and mortality. Note that the first particular distributed version of a classical quasispecies model was in fact introduced in [8, 9]. It was formulated as an integro-differential equation with respect to the sought-after time-varying probability density describing the relative frequencies of considered genotypes. A rigorous mathematical investigation of a class of such distributed models was provided in [29, 30, 31, 32, 25]. The dynamics in our distributed formulation is not restricted to the probability density constraint and takes growth and mortality characteristics into account.

Despite the mostly mathematical content of the current work, it has to be emphasized that the general methodology leading to our approach can be applied to a variety of different biological systems, including not only viruses but also bacteria, cancer cells, etc. Instead of descriptions of particular genotypes, more aggregated quantities may in principle be treated in quasispecies models. In particular, we keep in mind future applications of the proposed framework to developing mathematical models that involve directed therapy against various pathogens, such as infections or cancer [33, 34, 35], which will be considered elsewhere.

The rest of the paper is organized as follows. In Section 2, we construct a class of ODE-based open quasispecies models and establish its essential properties. Section 3 investigates the corresponding steady-state fitness maximization problem. Sections 2 and 3 also contain some related numerical simulation results. Section 4 is devoted to a class of continuous-time distributed open quasispecies models. A summarizing discussion of possible future developments is given in Section 5.

2 ODE-based open quasispecies models

2.1 Classical Eigen and Crow–Kimura quasispecies models

To make the text self-containing, we start with formulating the classical ODE-based quasispecies models.

The classical Eigen quasispecies model [3, 4] is formulated via the system of ordinary differential equations

[TABLE]

where the following notations are used:

• •

$t\geqslant 0$ is a time variable;

• •

$p(t)\>=\>(p_{1}(t),\,p_{2}(t),\,\ldots,\,p_{n}(t))^{\top}\>\in\>[0,+\infty)^{n\times 1}\>$ is the vector of the relative (normalized) time-varying frequencies of $n$ considered genotypes (or allelic effects) labeled by $1,2,\ldots,n$, $\>\sum_{i=1}^{n}p_{i}(t)\,=\,1\>$ for all $t\geqslant 0$;

• •

$K\,=\,\{k_{ij}\}_{i,j\,=\,1}^{n}\,\in\,[0,1]^{n\times n}\>$ is the mutation matrix, where $k_{ij}$ is the probability of the reproduction event for which an individual of genotype $j$ produces an individual of genotype $i$, so that

[TABLE]

• •

$M\>=\>\mathrm{diag}\,[m_{1},m_{2},\ldots,m_{n}]\>\in\>\mathbb{R}^{n\times n}\>$ is the diagonal matrix whose diagonal elements form the fitness landscape describing the relationship between the genotypes and reproductive success;

• •

$f[p(t)]\>=\>\sum_{i=1}^{n}m_{i}\,p_{i}(t)\>$ is the mean population fitness.

The Crow–Kimura quasispecies model [7] relies on the assumptions that the birth events are error-free and the mutations occur during the life time of the studied structures, i. e., the birth events and mutations are separated on the time scale. The corresponding dynamical system is written as

[TABLE]

where $\>\hat{M}\,=\,\{\mu_{ij}\}_{i,j\,=\,1}^{n}\,\in\,[0,+\infty)^{n\times n}$, $\>\mu_{ij}$ denotes the rate of mutation of genotype $j$ into genotype $i$,

[TABLE]

and other notations are as in (1). The permutation invariant Crow–Kimura model [36, 15] considers $N+1$ classes of genotypes and is obtained from (3) by replacing $n$ with $N+1$ and introducing the tridiagonal matrix

[TABLE]

By using the properties (2), (4) as well as the definitions of $M$ and $f$, one can easily verify that, for a solution of the system (1) or (3), the condition $\>\sum_{i=1}^{n}p_{i}(0)\,=\,1\>$ implies $\>\sum_{i=1}^{n}p_{i}(t)\,=\,1\>$ for all $t>0$.

The models (1) and (3) are stated for closed (isolated) systems of populations, such that their interactions with the environment and surrounding populations of other kinds are negligible with respect to the studied processes. The death rates are either absent or implicitly incorporated in the fitness landscape. Since the state variables are interpreted as relative frequencies and the state trajectories do not leave the standard simplex, growth and mortality are essentially not considered.

2.2 Open quasispecies models

From the perspective of mathematical modeling of therapy processes against infections or cancer [34, 33, 35], it is relevant to investigate how various quasispecies related to certain pathogens or diseased cells reproduce, mutate, die, and resist to targeted attacks from therapeutic agents. For that purpose, it is reasonable to build open quasispecies models whose state variables are interpreted not as relative frequencies, but as densities with respect to a certain external measure, so that their sum does not have to be invariant, and one can take growth and mortality properties into account. Before developing complicated models that involve the dynamics of a therapeutic agent and its influence on the dynamics of studied quasispecies, it is useful first to treat the simpler case when the death rates of the quasispecies are constant and the growth terms include only constant coefficients and a saturation factor depending on the total population size. This corresponds to the situation when the concentration of a therapeutic agent is assumed to be nearly constant within the observed time interval and can therefore be excluded from the state variables. Our work focuses namely on such kinds of open quasispecies models with growth and mortality (see also the discussion in the introduction).

Let us write both systems (1) and (3) in the common form

[TABLE]

with $\>G\,=\,\{g_{ij}\}_{i,j\,=\,1}^{n}\,\in\,\mathbb{R}^{n\times n}\>$ satisfying

[TABLE]

[TABLE]

One has $\,G=KM\,$ for (1) and $\,G=M+\hat{M}\,$ for (3).

Note that a solution of (6) with $\>\sum_{i=1}^{n}p_{i}(0)\,=\,1\>$ can be represented as

[TABLE]

where

[TABLE]

is the unique solution of

[TABLE]

[TABLE]

with an arbitrary constant $\>c\,\in\,\mathbb{R}\setminus\{0\},\>$ i. e., the components of $r$ may be understood as some size or density variables whose normalization leads to the relative frequencies constituting $p$ (a similar fact was originally noticed in [37]).

We hence proceed from the system (9) and incorporate growth and mortality characteristics in it. The following notations are adopted:

• •

$t\geqslant 0$ is a time variable;

• •

$u(t)\>=\>(u_{1}(t),\,u_{2}(t),\,\ldots,\,u_{n}(t))^{\top}\>\in\>[0,+\infty)^{n\times 1}\>$ is the vector of the densities of $n$ considered genotypes (or allelic effects) with respect to a certain external measure;

• •

$m\>=\>(m_{1},m_{2},\ldots,m_{n})^{\top}\>\in\>\mathbb{R}^{n\times 1}\>$ is the corresponding fitness landscape, $\>M\>=\>\mathrm{diag}\,\,m\>\in\>\mathbb{R}^{n\times n}$;

• •

$d\>=\>(d_{1},d_{2},\ldots,d_{n})^{\top}\>\in\>[0,+\infty)^{n\times 1}\>$ is the vector of the corresponding death rates, $\>D\>=\>\mathrm{diag}\,\,d\>\in\>[0,+\infty)^{n\times n}$;

• •

the growth and mutation properties are described by a matrix $G\in\mathbb{R}^{n\times n}$ satisfying (7), (8) as well as by a growth saturation term $\,\varphi\left(\sum_{i=1}^{n}u_{i}(t)\right)\,$ depending on the total population size $\,\sum_{i=1}^{n}u_{i}(t),\,$ where $\>\varphi\colon\,\mathbb{R}\to[0,+\infty)\>$ is an appropriate function.

Thus, we arrive at the autonomous system

[TABLE]

A similar way of introducing a growth saturation term was implemented in [38, 39] as applied to open replicator systems.

The origin $u=0$ is a trivial steady state of (10). Set the initial condition for (10) as

[TABLE]

By $\>F\>=\>(F_{1},F_{2},\ldots,F_{n})^{\top}\,\colon\>\mathbb{R}^{n\times 1}\to\mathbb{R}^{n\times 1},\>$ denote the function in the right-hand side of (10):

[TABLE]

If $\varphi$ is differentiable, the related Jacobian matrix takes the form

[TABLE]

Assumption 2.1**.**

$m$,* $M$, $d$, $D$, and $G$ are as given above in this subsection ((7) and (8) hold in particular*). $\varphi\colon\,\mathbb{R}\to[0,+\infty)\>$ is a nonnegative continuously differentiable function such that the functions $\>[0,+\infty)\,\ni\,s\>\longmapsto\>s\,\varphi(s)\>$ and $\>[0,+\infty)\,\ni\,s\>\longmapsto\>s\,\varphi^{\prime}(s)\>$ are bounded.

Theorem 2.2**.**

Let Assumption 2.1 hold. Then the nonnegative orthant $[0,+\infty)^{n\times 1}$ is positively invariant with respect to the system (10) (i. e., any state trajectory of (10) starting in this orthant stays there further up to the right end of the largest time interval of definition). Moreover,* for any initial state (11), there exists a unique solution of (10) that is defined for all $t\geqslant 0$.*

Theorem 2.2 is proved in Appendix.

Remark 2.3**.**

A particular admissible choice of the growth saturation function $\varphi$ is

[TABLE]

The case when $\>\varphi(s)\,=\,\alpha\,e^{-\gamma s}\>$ for all $s\in\mathbb{R}$ and $\alpha,\gamma$ are positive constants is trivially reduced to (14) by substituting $\>\varphi/\alpha,\,\alpha m,\,\alpha G\>$ with $\,\varphi,m,G,\,$ respectively. ∎

The next assumption is needed in particular to ensure the invertibility of the mortality matrix $D$ as well as the boundedness of the total population size function

[TABLE]

along the solution of (10), (11).

Assumption 2.4**.**

$d_{\min}\>=\>\min_{i\,=\,\overline{1,n}}\,d_{i}\>>\>0$.

Assumption 2.4 yields that

[TABLE]

Theorem 2.5**.**

Under Assumptions 2.1 and 2.4,* the total population size function (15) is bounded along the solution of (10), (11).*

Theorem 2.5 is proved in Appendix.

2.3 Steady-state analysis

In this subsection, we characterize a nontrivial steady state of the system (10) with growth and mortality under some new conditions in addition to Assumptions 2.1 and 2.4.

Assumption 2.6**.**

The function $\varphi$ is positive and strictly decreasing on $[0,+\infty)$,* and $\>l_{\varphi}\>=\>\lim_{s\,\to\,+\infty}\,\varphi(s)\>\geqslant\>0$.*

Remark 2.7**.**

Assumption 2.6 implies the existence of the strictly decreasing inverse function $\>\varphi^{-1}\,\colon\>(l_{\varphi},\varphi(0)]\,\to\,[0,+\infty)$. If $\varphi$ takes the form (14), one has

[TABLE]

and Assumption 2.6 is fulfilled. ∎

Recall also that the diagonal mortality matrix $D$ admits the diagonal inverse (16) in line with Assumption 2.4.

A nontrivial steady state $\>u^{*}\>=\>(u^{*}_{1},u^{*}_{2},\ldots,u^{*}_{n})^{\top}\>\in\>[0,+\infty)^{n\times 1}\>$ of the system (10) satisfies

[TABLE]

We therefore arrive at the problem of finding a positive real eigenvalue $\lambda^{*}$ of the matrix $D^{-1}G$ with a related eigenvector $\>u^{*}\>=\>(u^{*}_{1},u^{*}_{2},\ldots,u^{*}_{n})^{\top}\>\in\>[0,+\infty)^{n\times 1}\>$ such that

[TABLE]

For such $\lambda^{*}$ and $u^{*}$, the relations $\,Gu^{*}=\lambda^{*}Du^{*}\,$ and (7) lead to

[TABLE]

In the presence of mortality, it is reasonable to update the fitness definition for taking into account that an increase in the death rates negatively affects the population, as well as for maintaining a convenient representation of the steady-state fitness in terms of an appropriate eigenvalue similarly to the classical quasispecies models (for the latter, such a representation is discussed, e. g., in [25, §IV.3], [36, §2], [15, §1, §2], and [13, §1]).

Definition 2.8**.**

The population fitness for the modified quasispecies models described by the system (10) is given by

[TABLE]

Hence, the relation (19) means that

[TABLE]

One more assumption is required for investigating the eigenvalue problem for the matrix $D^{-1}G$. The definitions and properties of essentially nonnegative matrices, dominant eigenvalues, and irreducible matrices have to be recalled (see, e. g., [40, §I.7]).

Assumption 2.9**.**

The matrix $G$ is irreducible.

Remark 2.10**.**

The matrix $G$ is essentially nonnegative according to (8). Assumptions 2.4, 2.9 and the relation (16) imply that the matrix $D^{-1}G$ is also essentially nonnegative and irreducible (see [40, §I.7.4]). ∎

Remark 2.11**.**

For the matrix $\,G=KM\,$ in the Eigen model (1), Assumption 2.9 holds if, e. g., $\,m_{i}$, $i=\overline{1,n},\,$ and all elements of $K$ are positive (this sufficient condition can be relaxed in certain cases). For the matrix $\,G=M+\hat{M}\,$ in the Crow–Kimura model (3), Assumption 2.9 holds if and only if $\hat{M}$ is irreducible ($\hat{M}$ is a priori essentially nonnegative, since $M$ is diagonal and $\,G=M+\hat{M}\,$ fulfills (8)). In case of the permutation invariant Crow–Kimura model specified by (5), $\,\hat{M}=\mu Q\,$ is indeed irreducible, because $Q$ is irreducible. The latter property can be verified with the help of [40, §I.7.4] and [41, Fact 2 in §9.2] (the digraph based criterion of irreducibility for square matrices with nonnegative elements remains valid for essentially nonnegative matrices). ∎

Theorem 2.12**.**

Under Assumptions 2.1, 2.4,* and 2.9, the following properties hold*:**

• •

the matrix $D^{-1}G$ has a dominant eigenvalue $\lambda^{*}$, which is real and greater than the real part of any other eigenvalue of $D^{-1}G$;

• •

the eigenvalue $\lambda^{*}$ is simple (i. e., has algebraic multiplicity one)* and admits an eigenvector all of whose components are positive*;**

• •

there are no other eigenvalues of $D^{-1}G$ admitting eigenvectors all of whose components are nonnegative.

Proof.

It suffices to use Remark 2.10 and [40, Theorem I.7.10]. ∎

Theorem 2.13**.**

Let Assumptions 2.1, 2.4, 2.6,* and 2.9 hold. Then a nonzero steady state $\>u^{*}\>=\>(u^{*}_{1},u^{*}_{2},\ldots,u^{*}_{n})^{\top}\>\in\>[0,+\infty)^{n\times 1}\>$ of the system (10) exists if and only if the dominant eigenvalue $\lambda^{*}$ of the matrix $D^{-1}G$ satisfies*

[TABLE]

Moreover,* if (22) holds*,* this nontrivial steady state is uniquely determined as the eigenvector of $D^{-1}G$ corresponding to $\lambda^{*}$ and such that all of its components are positive and their sum equals $\,\varphi^{-1}(1/\lambda^{*})$*:**

[TABLE]

Proof.

It suffices to use Theorem 2.12 as well as the relations (17) and (18). ∎

Remark 2.14**.**

If $\varphi$ is selected according to (14), the relations (22) and (23) transform into

[TABLE]

(recall Remark 2.7). ∎

By using the well-known sufficient conditions for stability and instability from the first approximation (see, e. g., [42, Theorems 5.1 and 5.2 in Chapter 2]), it is easy to establish the following result.

Theorem 2.15**.**

Let Assumptions 2.1, 2.4, 2.6,* and 2.9 hold*,* and let $\>u^{*}\>=\>(u^{*}_{1},u^{*}_{2},\ldots,u^{*}_{n})^{\top}\>\in\>[0,+\infty)^{n\times 1}\>$ be a nonzero steady state of the system (10). Suppose also that the growth saturation function $\varphi$ is twice differentiable at the point $\>s^{*}\,=\,\sum_{i=1}^{n}u^{*}_{i}$. Consider the Jacobian matrix function $\mathrm{D}F$ determined by (13). At the steady state $u^{*}$, its elements can be represented as*

[TABLE]

If the real parts of all eigenvalues of $\,\mathrm{D}F(u^{*})$ are negative,* then the steady state $u^{*}$ is asymptotically stable. If at least one eigenvalue of $\mathrm{D}F(u^{*})$ has positive real part*,* then $u^{*}$ is unstable.*

Remark 2.16**.**

Recall the boundedness of the total population size along state trajectories as mentioned in Theorem 2.5. If a state trajectory is not attracted by the nontrivial steady state with positive components (see Theorem 2.13), the trivial steady state at the origin may be approached. The latter case means the extinction of the whole considered quasispecies population and is possible if all of the death rates are sufficiently large.

Note that the Jacobian matrix $\mathrm{D}F(0)$ at the origin equals $\>\varphi(0)\,G\,-\,D$. An interesting open problem is finding out sufficient conditions for the following dichotomic property which seems natural: either (i) the trivial steady state at the origin is asymptotically stable globally in the nonnegative orthant while the nonzero steady state $\,u^{*}\in[0,+\infty)^{n\times 1}\,$ does not exist, or (ii) the trivial steady state is unstable while the nontrivial steady state $u^{*}$ exists and is asymptotically stable globally in the nonnegative orthant with excluded origin. ∎

2.4 Numerical simulations

The numerical simulation results in this subsection demonstrate the possibility of two scenarios for our open models. Namely, with increasing a certain indicative parameter, the error threshold is observed when the population does not extinct, or, alternatively, the extinction already occurs prior to the error threshold. In the second case, the error threshold is understood nominally.

For simplicity, we adopt the exponential form (14) of the growth saturation function $\varphi$ and consider the dynamical system (10) with the permutation invariant Crow–Kimura formalism such that $n=N+1$, $\,G=M+\hat{M},\,$ and $\hat{M}$ is given by (5). The selected form of $\varphi$ allows for using the simplifications mentioned in Remarks 2.7 and 2.14.

We take

[TABLE]

while the mutation rate parameter $\mu$ (see (5)) and the death rate $\Delta$ for the genotypes of the zero class are not fixed, so that the key quantities $\lambda^{*}$ and $u^{*}$ become functions of $\mu$ and $\Delta$.

Fig. 1 illustrates the situation when, with the increase of the mutation rate $\mu$, the population clearly goes through the error threshold with no extinction. It is also demonstrated that, as the death rate $\Delta$ of the zero class grows, the error threshold occurs for lower $\mu$. The graphs of the maximum over the real parts of the eigenvalues of $\mathrm{D}F(u^{*})$ yield asymptotic stability of $u^{*}$ for all observed values of $\mu$ except for the critical value of $\mu$ specifying the error threshold. However, it is natural to expect the stability for this critical $\mu$ as well. Besides, extensive numerical simulations allow for conjecturing that the existence of the nontrivial steady state $u^{*}$ would imply its asymptotic stability globally in the nonnegative orthant with excluded origin (recall Remark 2.16).

Fig. 2 shows that the death rate $\Delta$ can also be considered as a bifurcation parameter whose changes eventually lead to the error threshold.

We now change the death rate vector to

[TABLE]

(i. e., the death rates of all classes except for the zero class are increased from $0.5$ to $2$), while all other parameters in (26) remain the same. Fig. 3 illustrates that the error threshold nominally appears already after the extinction threshold $\lambda^{*}=1$ and hence has almost no effect on the dynamics. The extinction takes place when $0<\lambda^{*}\leqslant 1$; nonpositive values of $\lambda^{*}$ are not observed. Note that, according to Theorem 2.13 and Remark 2.14, the nontrivial steady state $u^{*}$ exists if and only if $\lambda^{*}>1$. For $0<\lambda^{*}\leqslant 1$, the right-hand side of the normalization condition $\>\sum_{i=1}^{n}u^{*}_{i}\,=\,\ln(\lambda^{*})/\gamma\>$ (see (24)) becomes negative, and the nontrivial steady state does not exist, so we understand $u^{*}$ nominally.

The presented numerical simulation results indicate that even the simplified permutation invariant Crow–Kimura model in our open setting possesses a rich dynamical behavior with different bifurcation scenarios.

3 Steady-state fitness maximization

3.1 Theoretical analysis

Similarly to the discussion in [43, Introduction], we adopt the hypotheses that evolutionary adaptation of the fitness landscape of a quasispecies population is significantly slower than the internal dynamics of the corresponding state variables (in our case, this dynamics is governed by the system (10)), and that changes in the fitness landscape are aimed at eventually maximizing the population fitness (specified by Definition 2.8 in our case) under some constraints. Based on these hypotheses as well as on the relations (20), (21) and Theorem 2.13, we arrive at the problem of maximizing the population fitness

[TABLE]

at the nonzero steady state $\>u^{*}[m]\>=\>(u^{*}_{1}[m],\,u^{*}_{2}[m],\,\ldots,\,u^{*}_{n}[m])^{\top}\>\in\>[0,+\infty)^{n\times 1}\>$ of (10) over fitness landscapes $\>m\>=\>(m_{1},m_{2},\ldots,m_{n})^{\top}\>\in\>\mathbb{R}^{n\times 1}\>$ subject to the constraints

[TABLE]

[TABLE]

where $\Pi$ is a bounded subset of $\mathbb{R}^{n\times 1}$, and the death rates $\,d_{i}$, $i=\overline{1,n},\,$ are fixed for simplicity.

As was mentioned in the introduction, we try to keep the connection with the classical quasispecies models as close as reasonably possible. For our open quasispecies models, it is therefore useful to consider first the simplified case when mutations affect the fitness landscape (and, consequently, the growth rates), but not the death rates. Modeling the influence of mutations also on the death rates is a possible subject of future research.

Besides, it is useful to note the paper [44] that contains a relevant discussion on the applicability of the fitness maximization principle, though as applied to classical replicator systems.

As was discussed in Remark 2.16, the nontrivial steady state may be unstable for some admissible fitness landscapes (when the growth terms are overriden by the mortality terms). However, one would expect that, if any optimal steady state is unstable, all other admissible steady states will also be unstable, although a rigorous derivation of sufficient conditions for this property remains an open problem. Note also that, in case $l_{\varphi}=0$, the constraint (29) is simplified to

[TABLE]

For example, if $\varphi$ is given by (14), then (29) transforms into $\,\lambda^{*}[m]>1\,$ (see Remark 2.7).

Thus, the problem of maximizing the dominant eigenvalue $\lambda^{*}[m]$ of the matrix $\,D^{-1}\cdot G[m]\,$ (see Theorem 2.13), i. e.,

[TABLE]

plays a central role.

Remark 3.1**.**

For the classical quasispecies models described by the system (6), one has to consider the dominant eigenvalue of the matrix $G$, and, therefore, the subsequent investigation regarding the optimization problem (31) will also be valid in the classical case after the formal replacement of $D$ with the identity matrix of size $n\times n$. ∎

Remark 3.2**.**

If $\lambda^{*}[m]$ is defined for all $m\in\Pi$, then this is a continuous function of $m\in\Pi$ due to the well-known fact that the eigenvalues of a real or complex square matrix continuously depend on its entries (see, e. g., [45, Theorem 2.11]). ∎

An important case is when (31) becomes a convex optimization problem.

Assumption 3.3**.**

$\Pi$* is a convex compact subset of $\,\mathbb{R}^{n\times 1}$.*

Theorem 3.4**.**

Let Assumptions 2.1, 2.4, 2.6, 2.9,* and 3.3 hold*,* and let the matrix $G$ have the Crow–Kimura form $\,G=M+\hat{M},\,$ where $\,M\,=\,\mathrm{diag}\,m,\,$ while $\hat{M}$ does not depend on $m$. Then the function*

[TABLE]

is convex.

Proof.

It suffices to use the result of [46] saying that the dominant eigenvalue of an essentially nonnegative matrix is convex if considered as a function of the diagonal elements of this matrix. ∎

Remark 3.5**.**

If $G$ has the Eigen form $\,G=KM,\,$ where $\,M\,=\,\mathrm{diag}\,m,\,$ while $K$ satisfies (2) and does not depend on $m$, then the result of [46] cannot be applied directly. In this case, obtaining sufficient conditions for the convexity of the dominant eigenvalue $\lambda^{*}[m]$ of $\,D^{-1}\cdot G[m]\,$ with respect to $m$ remains an open problem. ∎

We now assume that the considered system satisfies this convexity property on the fitness landscape constraint set $\Pi$.

Assumption 3.6**.**

The matrix $\,D^{-1}\cdot G[m]\,$ admits the dominant eigenvalue $\,\lambda^{*}[m]\in\mathbb{R}\,$ for any $m\in\Pi$,* and the function*

[TABLE]

is convex.

Recall that a point $x$ of a convex set $C\subseteq\mathbb{R}^{n}$ is called an extreme point of $C$ if and only if there is no way to express $x$ as a convex combination $\>\alpha y\,+\,(1-\alpha)z\>$ with $y\in C$, $z\in C$, and $\alpha\in(0,1)$, except by putting $y=z=x$ (see, e. g., [47, §18]).

Theorem 3.7**.**

Let Assumptions 2.1, 2.4, 2.6, 2.9, 3.3,* and 3.6 hold. Then the function (32) is continuous on $\Pi$ and attains its maximum over $\Pi$, and any related maximizer is an extreme point of $\Pi$.*

Proof.

It suffices to take Remark 3.2, Assumption 3.3, and [47, Corollary 32.3.1] into account. ∎

Corollary 3.8**.**

If the conditions of Theorem 3.7 hold and,* moreover*,* $\Pi$ is a compact convex polytope*,* then any maximizer in the problem (31) is a vertex of $\Pi$.*

3.2 Numerical simulations

The numerical simulation results in this subsection are purely illustrative and demonstrate a typical behavior of the maximal steady-state population fitness with the increase of a mutation rate or a death rate.

We adopt the exponential form (14) of the growth saturation function $\varphi$ and consider the dynamical system (10) with the following symmetric Crow–Kimura formalism [48, 44]:

• •

$n=2^{N}$ genotypes are considered, they are labelled by $\,0,1,\ldots,n-1\,$ and associated to the $N$-dimensional binary representations (sequences) of the corresponding numbers;

• •

$G\,=\,M+\hat{M}$, $\>M\,=\,\mathrm{diag}\,m$, $\>\hat{M}\,=\,\{\hat{m}_{ij}\}_{i,j\,=\,1}^{n}$,

[TABLE]

where $\mu>0$ is a mutation rate parameter and $H_{ij}$ denotes the Hamming distance between sequences $i$ and $j$ ($H_{ij}=0$ only if $i=j$).

The irreducibility of such a mutation matrix $\hat{M}$ can be verified via the same reasonings as in the end of Remark 2.11 for the permutation invariant Crow–Kimura model. Hence, Theorems 3.4, 3.7 and Corollary 3.8 can be applied, and the steady-state fitness maximization problem for fitness landscapes constrained by a compact convex polytope $\Pi$ reduces to sorting through the vertices of $\Pi$.

Here we do not use the permutation invariant Crow–Kimura formalism (see (5)) because of its high heterogeneity such that many of the original $2^{N}$ binary sequences are concentrated in the middle of the ordered group of the $N+1$ equivalence classes. Thus, for attenuating the related heterogeneity in the steady-state fitness profile, one typically has to select a rather specific structure of the death rate vector.

For simplicity, we consider the low-dimensional case

[TABLE]

Moreover, we take

[TABLE]

as in (26). The components of the death rate vector $d$ are independently and randomly generated from the uniform distribution on $(0,1)$, and they are also arranged in ascending order:

[TABLE]

The constraint set $\Pi$ is chosen as the following simplex with center at the point $\>(0.5,0.5,\ldots,0.5)^{\top}\,\in\,\mathbb{R}^{n\times 1}$:

[TABLE]

There are $n=8$ vertices which we label by $\>1,2,\ldots,n$: for every $\>i\,\in\,\{1,2,\ldots,n\},\>$ the $i$-th coordinate of vertex $i$ is $2.5$, while all other coordinates of vertex $i$ are equal to $0.5$.

Fig. 4 indicates how the steady-state fitnesses at vertices $1,2,3$ depend on the mutation rate parameter $\mu$. The maximal steady-state fitness in the problem (31) is achieved at vertex 1. This is not surprising, since the components of the death rate vector (36) are written in ascending order. Furthermore, when $\mu$ is sufficiently large, the difference between the steady-state fitnesses at the vertices of $\Pi$ (not only 1–3, but also 4–8) becomes negligible.

Next, we fix $\mu$ and let the first component of the death rate vector be a variable $\Delta$:

[TABLE]

Fig. 5 illustrates how the steady-state fitnesses at vertices 1 and 2 depend on the death rate $d_{1}=\Delta$ of genotype 0 (our notation is such that, for every $\>i\,\in\,\{1,2,\ldots,n\},\>$ the $i$-th component of $d$ is the death rate of genotype $i-1$). The maximal steady-state fitness in the problem (31) is achieved at vertex 1 if $\,\Delta\leqslant d_{2}=0.158373\,$ and at vertex 2 if $\Delta\geqslant d_{2}$.

4 Continuous-time distributed open quasispecies models

4.1 Constructing the dynamical equation

Since the actual number of genotypes or alleles involved in a studied process can be extremely large (even though a significant fraction of them may have small effects), distributed (continuum-of-alleles) quasispecies models are of particular interest (see, e. g., [25, Introduction to Chapter IV]). For the classical Eigen and Crow–Kimura models, distributed extensions in the form of an integro-differential equation with respect to a sought-after time-varying probability density were rigorously investigated in [29, 30, 31, 32] and [25, Chapter IV], while the first particular distributed model of this kind was proposed in [8, 9]. Bürger’s formalism treats the distributed Eigen and Crow–Kimura settings together and allows to establish sufficient conditions for the existence and uniqueness of an equilibrium probability density with its positivity and asymptotic stability (see [30, §3, §4] and [25, §IV.3]).

Typical (although not the most general) integro-differential equations describing the Eigen and Crow–Kimura dynamics can be written as

[TABLE]

and

[TABLE]

respectively, where the following notations are used:

• •

$t\geqslant 0$ is a time variable;

• •

the considered genotypes (or allelic effects) are represented as points of a finite-dimensional region $\Omega$;

• •

$p\,\colon\>\Omega\times[0,+\infty)\>\to\>[0,+\infty)\>$ is the function describing the relative frequencies of the genotypes and normalized so that

[TABLE]

• •

$\Omega\,\ni\,x\>\longmapsto\>K(x,y)\,\in\,[0,+\infty)\>$ is the probability density related to reproduction events such that an individual of genotype $y\in\Omega$ produces an individual of genotype $x\in\Omega$, i. e., $\>K\colon\,\Omega^{2}\to[0,+\infty)\>$ is a distributed analog of the mutation matrix, and

[TABLE]

• •

$m\colon\,\Omega\to\mathbb{R}\>$ is the distributed fitness landscape;

• •

$f[p(\cdot,t)]\>=\>\int_{\Omega}m(x)\,p(x,t)\,\mathrm{d}x\>$ is the mean population fitness;

• •

$\mu\colon\,\Omega\,\to\,[0,+\infty)\>$ is the mutation rate function.

For simplicity, we impose some more specific conditions than in [30, §3, §4] and [25, §IV.3]. In particular, it is reasonable to consider a bounded region $\Omega$ from a computational perspective.

Assumption 4.1**.**

$\Omega$* is a bounded open domain in $\mathbb{R}^{\varkappa}$ or the closure of such a domain*,* $\varkappa\in\mathbb{N}$, the functions $\>K\colon\,\Omega^{2}\to[0,+\infty)$, $\>m\colon\,\Omega\to\mathbb{R},\>$ and $\>\mu\colon\,\Omega\,\to\,[0,+\infty)\>$ are bounded*,* and (42) holds.*

One can consider (39) and (40) as differential equations whose state space is the Banach space

[TABLE]

of all Lebesgue measurable functions $\>v\colon\,\Omega\to\mathbb{R}\>$ for which the Lebesgue integral $\>\int_{\Omega}|v(x)|\,\mathrm{d}x\>$ exists and is finite (the theory of differential equations in Banach spaces is introduced, e. g., in [49]).

The relations $\,p(x,t)\geqslant 0\,$ and (41) are preserved along the solutions of (39) and (40).

Both equations (39) and (40) can be written in the common form

[TABLE]

with the functions $\>G\colon\,\Omega^{2}\to\mathbb{R}\>$ and $\>b\colon\,\Omega\to\mathbb{R}\>$ given by

[TABLE]

or

[TABLE]

Note that a solution of (44) with $\>\int_{\Omega}p(x,0)\,\mathrm{d}x\,=\,1\>$ can be represented as

[TABLE]

where

[TABLE]

is the unique solution of

[TABLE]

[TABLE]

with an arbitrary constant $\>c\,\in\,\mathbb{R}\setminus\{0\},\>$ i. e., the values of $r$ may be interpreted as some quantities whose normalization leads to the values of the time-varying probability density $p$.

We therefore proceed from the system (47) and incorporate growth and mortality characteristics in it (recall similar considerations in Subsection 2.2). Let $\>u\,\colon\>\Omega\times[0,+\infty)\>\to\>[0,+\infty)\>$ be a function specifying the sought-after dynamical quantities associated with the considered genotypes. Introduce also the death rate function $\>d\colon\,\Omega\,\to\,[0,+\infty)\>$ and the growth saturation term $\>\varphi\left(\int_{\Omega}u(x,t)\,\mathrm{d}x\right)$.

Assumption 4.2**.**

$d\colon\,\Omega\,\to\,[0,+\infty)\>$* is bounded*,* $\>\varphi\colon\,\mathbb{R}\to[0,+\infty)\>$ fulfills the same properties as mentioned in Assumption 2.1, $\>G\colon\,\Omega^{2}\to\mathbb{R}\>$ and $\>b\colon\,\Omega\to\mathbb{R}\>$ are defined by either (45) or (46), and*,* moreover*,**

[TABLE]

Thus, we arrive at the integro-differential equation

[TABLE]

Consider it as a differential equation with the state space (43), and set the initial condition as

[TABLE]

where

[TABLE]

The zero element in $V$ is a trivial steady state of (48).

For any real Banach spaces $W_{1}$ and $W_{2}$, let $\mathcal{L}(W_{1},W_{2})$ denote the set of all bounded linear operators acting from $W_{1}$ to $W_{2}$.

It is convenient to rewrite the right-hand side of (48) by means of the following operator notation:

[TABLE]

Theorem 4.3**.**

Under Assumptions 4.1 and 4.2,* the Fréchet derivative (Jacobian operator) $\>\mathrm{D}\mathcal{F}\,\colon\>V\,\to\,\mathcal{L}(V,V)\>$ of the right-hand side operator $\mathcal{F}$ is represented as*

[TABLE]

where $\,\mathcal{E}_{v}\,\in\,\mathcal{L}(V,V)\,$ is defined by

[TABLE]

for any $v\in V$.

Theorem 4.4**.**

Let Assumptions 4.1 and 4.2 hold. Then the subset $V_{+}$ of the state space $V$ is positively invariant with respect to the dynamical system (48). Moreover,* for any initial state (49), there exists a unique solution of (48) that is defined for all $t\geqslant 0$.*

Theorems 4.3 and 4.4 are proved in Appendix.

The next condition is imposed in particular to guarantee the boundedness of the total population size function

[TABLE]

along the solution of (48), (49).

Assumption 4.5**.**

There exists a constant $d_{\mathrm{low}}>0$ such that $\,d(x)\geqslant d_{\mathrm{low}}\,$ for all $x\in\Omega$.

Theorem 4.6**.**

Under Assumptions 4.1, 4.2,* and 4.5, the total population size function (54) is bounded along the solution of (48), (49).*

Theorem 4.6 can be proved similarly to Theorem 2.5. Compared to the proof of the latter (given in Appendix), sums over indices $\,1,2,\ldots,n\,$ should now be replaced with integrals over $\Omega$, while $m_{\max}$ and $d_{\min}$ should be replaced with $\>\mathrm{ess}\,\sup_{x\,\in\,\Omega}\,m(x)\>$ and $d_{\mathrm{low}}$, respectively.

4.2 Steady-state analysis

For investigating a nontrivial steady state of the integro-differential equation (48), additional conditions have to be adopted.

Assumption 4.7**.**

The following properties hold:**

the function $\varphi$ satisfies the same conditions as mentioned in Assumption 2.6;

• 2)

for every measurable set $\Omega_{1}\subset\Omega$ such that the Lebesgue measures of both $\Omega_{1}$ and $\Omega\setminus\Omega_{1}$ are positive,* one has*

[TABLE]

• 3)

in the Eigen case (45), there exists a subset $\hat{\Omega}\subseteq\Omega$ of positive Lebesgue measure and a constant $\hat{c}>0$ such that

[TABLE]

• 4)

in the Crow–Kimura case (46), there exists a subset $\hat{\Omega}\subseteq\Omega$ of positive Lebesgue measure and a constant $\hat{c}>0$ such that (55) holds together with

[TABLE]

and

[TABLE]

where divergence of the integral is allowed.

Remark 4.8**.**

Let Assumptions 4.1, 4.2, and 4.5 hold. Then Item 2 of Assumption 4.7 holds if, e. g., $G$ is positive on $\Omega^{2}$. Item 2 means that the compact operator $\,\mathcal{G}_{1}\,\in\,\mathcal{L}(V,V)\,$ given by

[TABLE]

is irreducible (see [30, Proposition 3.1]). Regarding Item 4, the integral in (56) diverges if, e. g., there exist a point $\hat{x}\in\hat{\Omega}$ and a constant $c_{1}>0$ such that

[TABLE]

∎

A nontrivial steady state $u^{*}\in V_{+}$ of the integro-differential equation (48) satisfies

[TABLE]

where the operator notations (57) and

[TABLE]

are used. We hence arrive at the problem of finding a positive real eigenvalue $\lambda^{*}$ of the operator $\>\mathcal{G}_{1}+\mathcal{B}_{1}\,=\,\mathcal{D}^{-1}(\mathcal{G}+\mathcal{B})\>$ with a related eigenfunction $u^{*}\in V_{+}$ such that

[TABLE]

(recall similar reasonings in Subsection 2.3 with the formulas (17) and (18)). For such $\lambda^{*}$ and $u^{*}$, the relations $\>\mathcal{G}[u^{*}]\,=\,\lambda^{*}\,\mathcal{D}[u^{*}]\>$ and (42), (45), (46), (51) lead to

[TABLE]

(recall also (19)).

Similarly to Definition 2.8, we specify the fitness function as follows.

Definition 4.9**.**

The population fitness for the modified distributed quasispecies models described by the integro-differential equation (48) is given by

[TABLE]

The relation (61) therefore means that

[TABLE]

Theorem 4.10**.**

Under Assumptions 4.1, 4.2, 4.5,* and 4.7, the following properties hold*:**

• •

the operator $\>\mathcal{G}_{1}+\mathcal{B}_{1}\,=\,\mathcal{D}^{-1}(\mathcal{G}+\mathcal{B})\>$ has a dominant eigenvalue $\lambda^{*}$, which is real and greater than the real part of any other spectral point of $\mathcal{G}_{1}+\mathcal{B}_{1}$;

• •

the eigenvalue $\lambda^{*}$ is simple and admits a positive eigenfunction;**

• •

there are no other eigenvalues of $\mathcal{G}_{1}+\mathcal{B}_{1}$ admitting nonnegative eigenfunctions.

Theorem 4.10 is proved in Appendix.

Theorem 4.11**.**

Let Assumptions 4.1, 4.2, 4.5,* and 4.7 hold. Then a steady state $u^{*}\in V_{+}$ of the integro-differential equation (48) with $\>s^{*}\>=\>\int_{\Omega}u^{*}(x)\,\mathrm{d}x\>>\>0\>$ exists if and only if the dominant eigenvalue $\lambda^{*}$ of the operator $\mathcal{G}_{1}+\mathcal{B}_{1}$ satisfies*

[TABLE]

Moreover,* if (64) holds*,* this steady state is uniquely determined as the positive eigenfunction of $\mathcal{G}_{1}+\mathcal{B}_{1}$ corresponding to $\lambda^{*}$ and normalized so that*

[TABLE]

Proof.

It suffices to use Theorem 4.10 as well as the relations (58) and (60). ∎

Note the similarity between Theorems 4.10, 4.11 and 2.12, 2.13, respectively.

Remark 4.12**.**

If $\varphi$ is given by (14), the relations (64) and (65) transform into

[TABLE]

(recall Remark 2.14). ∎

Next, let us provide sufficient conditions for stability and instability of the nontrivial steady state $u^{*}$ of (48). By using the general results mentioned in [49, §VII.2.4], one obtains the following theorem.

Theorem 4.13**.**

Let Assumptions 4.1, 4.2, 4.5,* and 4.7 hold*,* and let $u^{*}\in V_{+}$ satisfy $\>s^{*}\>=\>\int_{\Omega}u^{*}(x)\,\mathrm{d}x\>>\>0\>$ and be a steady state of the integro-differential equation (48). Suppose also that the growth saturation function $\varphi$ is twice differentiable at the point $s^{*}$. Consider the Jacobian operator $\mathrm{D}\mathcal{F}$ determined by (52) and (53). The auxiliary operator $\mathcal{E}_{u^{*}}$ for the steady state $u^{*}$ is simplified to*

[TABLE]

If the real parts of all spectral points of $\mathrm{D}\mathcal{F}[u^{*}]$ are negative,* then the steady state $u^{*}$ is asymptotically stable. If at least one spectral point of $\mathrm{D}\mathcal{F}[u^{*}]$ has positive real part*,* then $u^{*}$ is unstable.*

Remark 4.14**.**

If $\varphi$ is selected in line with (14), the relation (67) transforms into

[TABLE]

∎

Remark 4.15**.**

Recall the boundedness of the total population size along state trajectories as mentioned in Theorem 4.6. If a state trajectory is not attracted by the positive steady state, the zero steady state may be approached. The latter case means the extinction of the whole considered quasispecies population and is possible when $\>\mathrm{ess}\,\inf_{x\,\in\,\Omega}\,d(x)\>$ is sufficiently large (see also Remark 2.16). ∎

The challenging practical problems of obtaining the dominant eigenvalue $\lambda^{*}$ and the positive eigenfunction $u^{*}$ of the operator $\mathcal{G}_{1}+\mathcal{B}_{1}$ and verifying stability or instability of the steady state $u^{*}$ via Theorem 4.13 (so that one has to characterize the spectrum of $\mathrm{D}\mathcal{F}[u^{*}]$) are possible subjects of future research. Their principal issue is that eigenvalue problems for noncompact and non-self-adjoint compact operators are involved. In general, the compact integral operator $\mathcal{G}_{1}$ is non-self-adjoint (due to its nonsymmetric kernel), the operator $\mathcal{G}_{1}+\mathcal{B}_{1}$ is compact in the Eigen case (45) (when $\mathcal{B}_{1}$ vanishes) but noncompact in the Crow–Kimura case (46), and $\mathrm{D}\mathcal{F}[u^{*}]$ is noncompact in both of these cases (since $\mathcal{D}$ is noncompact). After an efficient approach to treating such problems is developed (at least for a particular nontrivial subclass of the distributed quasispecies models), it will be reasonable to investigate the infinite-dimensional problems of maximizing the steady-state fitness over fitness landscape functions that satisfy appropriate constraints.

5 Conclusion

In this work, we developed a general approach to the construction of open quasispecies models incorporating growth and mortality characteristics. We proposed open modifications of the Eigen and Crow–Kimura quasispecies models in both ODE-based and distributed (continuum-of-alleles) settings. The distributed formulations built on integro-differential equations and were motivated by the complexity of numerical analysis of quasispecies systems with large numbers of ODEs, as well as by the fact that the actual number of genotypes or alleles involved in a studied process could indeed be extremely large. Essential properties of the open quasispecies models, regarding in particular steady states, were investigated.

We also explained a motivation for steady-state fitness maximization problems and studied them in case of ODE-based quasispecies dynamics. It was in particular established that, for the ODE-based Crow–Kimura models (both open and closed), such a problem leads to convex optimization and allows for an efficient numerical implementation. For the Eigen models, it remains an open question whether a similar reduction to convex optimization can in general be carried out.

Another open problem for our models is verification of the following natural conjecture under the already adopted and possibly some additional assumptions: either (i) the trivial zero steady state is asymptotically stable globally in the nonnegative subspace while the nontrivial nonnegative steady state does not exist, or (ii) the trivial steady state is unstable while the nontrivial steady state exists and is asymptotically stable globally in the nonnegative subspace with excluded zero.

As was demonstrated in our numerical simulation results, the ODE-based open quasispecies models enable at least the following two nontrivial bifurcation scenarios with the increase of a specific parameter (such as a mutation rate or a death rate): (i) the error threshold is observed when the population does not extinct, or (ii) the extinction already takes place prior to the nominally interpreted error threshold.

Our mathematical constructions for the distributed open quasispecies models serve just as the first step in this research direction. Further developments require the design of an efficient practical approach to treating the related eigenvalue problems and verifying stability or instability of the nontrivial nonnegative steady states. A major difficulty in these problems is that one has to deal with noncompact and non-self-adjoint compact operators. After the corresponding challenges are overcome at least for a particular reasonable subclass of the distributed open quasispecies models, it will be relevant to consider the infinite-dimensional problem of steady-state fitness maximization.

As was also noted above, a promising research area where our framework may eventually be used is modeling the dynamics of various quasispecies related to certain pathogens or diseased cells under targeted attacks from therapeutic agents [33, 34, 35].