1. Introduction

Epidemiological studies concern with models of the occurrence of diseases in space and time. They are carried out relying on mathematical models and their analysis' tools; mathematic aids in inferring disease causes, predicting the future course of an outbreak and planning the most appropriate control measures. The most numerous epidemiological mathematical models are in the field of infectious diseases [11,19,24]. Under different assumptions, dynamical systems have been designed to predict the number of infected people over times, as well as to quantify the effect of diseases' spread control measures such as pharmacological therapies or social actions. [2,22,41].

When referring to genetic diseases, epidemiological studies examine the role of genetic factors and their interaction with environmental features in the occurrence of human genetic diseases i.e., an illness caused by abnormalities within the DNA sequence of human genes. Their general aim is to understand the role of genetic factors in the etiology of diseases in human populations with the ultimate goal of disease control and prevention [25].

Mathematical models and analyses have successfully been applied in human genetics epidemiology (see [25,26,27] and reference therein). Dynamic models have been developed to describe gene distributions in evolving populations, as well as changes in a single genome or biochemical process regulated by genes [16,43,17]. Mathematical models can help to provide non intuitive insights into disease spread within a population as well as to predict the influence of de novo gene mutations on the incidence of these diseases or quantify the effectiveness of treatments [7,3,44].

We studied the epidemiology of a class of genetic disorders, namely the X-linked recessive diseases; these include the serious diseases hemophilia A, Duchenne/Becker muscular dystrophy, and Lesch-Nyhan syndrome as well as common and less serious conditions such as male pattern baldness and red-green colour blindness. Although genetic diseases are very rare, it has been estimated that millions of people are affected worldwide [21]. The incidence of an X-linked recessive disease depends on the disorder's severity: it ranges from 1 in 3000 newborn males for Duchanne/Becker muscular dystrophy to 1 in 20 for the red and green color blindness.

Despite their relevance, only a few mathematical models have been specifically developed to study the dynamics and spread of X-linked recessive diseases in a population [43,14]. Related studies analyze the inheritance mechanism of any gene-not necessarily responsible for a genetic disease-placed on the X chromosome [29,28]; they belong to the field of population genetics. In these works genotypes frequencies-i.e., the frequency or proportion of genotypes in a population-are often chosen as model's variables. Under the hypothesis of infinite population and starting from a genotypes' distribution, the genotypes' proportions in the next generation are evaluated according to the inheritance mechanism and to the effects of selection or mutation. The average fitness (see [30] pag. 385-387) is frequently studied as a suitable Lyapunov function candidate to analyze stability properties of model's equilibrium points. Even in this generic scenario seldom contributions examine the combined effects of selection and mutation on population's dynamics and equilibrium (see [31,39] and reference therein). Furthermore findings of these researches can not be applied to epidemiological studies. In fact the ultimate goal of genetic epidemiology is to predict the number of individuals carrying the disease responsible gene; this number can not be inferred from genotypes' frequency distribution when assuming infinite population size.

We developed a discrete time mathematical model aiming at describing the spread of X-linked recessive diseases in the population as well as quantifying the influence on diseases distribution of related factors such as individual's fitness and sporadic genetic mutation.

In our preliminary studies [14,15] we developed an X-linked recessive diseases model which did not take into account the class of affected women. In [14] we made the simplifying hypothesis that only sons born from healthy couples could be born affected owing to spontaneous mutations. Fitness rates were assumed identical for all individuals, regardless of their health conditions. In that simplified scenario we proved that the population asymptotically distributes among classes in proportions depending on the spontaneous mutation rate values. In [15] we studied the influence of fitness factors on the same model framework and we inferred general system properties assuming that either selection or mutation apply.

Encouraged by previous results and leveraging on the same framework, we elaborated the model we are here presenting which includes the class of affected women leading to a five state system. We also took into account genetic mutations possibly occurring to any couples' offspring and fitness factors varying according to individual's health conditions. These features enable to reproduce the inheritance mechanism and the spread in a population of any X-linked disease.

Although the mathematical model we developed is nonlinear, it is suitable to certain analyses using classical nonlinear methods to gain information about system behavior, equilibrium existence and their convergence properties. We remark that in genetic epidemiology convergence means that a population initialized at some point (i.e., with an assigned initial allocation among healthy, affected and carrier people) will move, through a time sequence of generations, to a (possibly) different population distribution [10]. Similar analysis have been conducted in general biological systems: in [4] a Lyapunov method is proposed for finding the invariant sets of non-negative dynamical systems modeling overpopulation species in ecology, while in [8] equivalent results are obtained using Jacobian analysis.

The organization of the paper is as follow. In Section 2 we introduce a background on X-linked recessive diseases, details on spontaneous genetic mutation and individual's fitness factors. After giving details on the assumptions and on the inheritance mechanisms of X-linked diseases we propose a mathematical model in Section 3 and its equilibrium points are analyzed in Section 4. Simulation results and a sensitivity analysis are given respectively in Section 5 and 6. Finally, some discussions and conclusions are presented in Section 7 and 8.

---

2. The X-linked recessive inheritance mechanism

Clinical expression and inheritance patterns of X-linked recessive diseases are related to individual gender. While females possess two X chromosomes, males have one X chromosome and one Y chromosome; thus females possess two copies of each X chromosome gene, whereas males only have a single copy of each X and Y chromosome genes [16].

X-linked diseases can be either recessive or dominant [18]. In X-linked recessive diseases females harbor one copy of each the disease responsible and the normal gene. As the effects of the abnormal gene are counteracted by those of the normal gene located on the alternative X chromosome, females usually show no signs of disease. For that reason, they are called carriers, and can transmit either the normal or the abnormal gene to their progeny. Thus, to become affected by an X-linked recessive disease, females have to harbor a copy of the disease-related gene on both X chromosomes. In contrast, males harboring an abnormal gene on their unique X chromosome are usually affected, even if the disease is recessive. Indeed, the Y chromosome lacks all X chromosome genes, and is therefore unable to exert any compensatory effect.

An X-linked recessive disease may be inherited according to the following rules:

• Affected males never spread the disease to their sons, as no male-to-male transmission of the X chromosome occurs (see Figure 1(a)).

• Affected males pass the abnormal X chromosome to all of their daughters, who are described as obligate carriers (see Figure 1(a)).

• On average, female carriers pass a defective X chromosome to half of their sons (who will born affected) and half of their daughters (who become carriers). The remaining half of their siblings inherit a normal copy of the chromosome (see Figure 1(b)).

• Affected females are the rare result of an affected male and a carrier or affected female mating. They are more frequent in less serious X-linked recessive diseases that do not comprise reproduction capacities.

We refer the reader to [18] for details on clinical expression and inheritance patterns of X-linked recessive diseases.

Other than the result of the described, well characterized patterns of inheritance, the current spread of genetic disorders within a given population is influenced by additional factors, including sporadic mutations and individual's fitness. Both are "driving forces" in the transmission mechanism of X-linked recessive diseases and they will be described in the next section; other factors such as prenatal diagnosis, population migrations and non random mating slightly affect diseases' spread [31,35].

---

2.1. Spontaneous genetic mutation

A genetic disease that occurs when neither parent is affected or carrier of any genetic defect is called sporadic mutation or de novo gene mutations. These cases arise via random genetic mutations within the DNA sequence; the mutation can occur in the germ-line cell population-i.e. in eggs and sperm cells-in subjects without any prior genetic defect and it can be transmitted down to one of the offspring [5]. The genetic mutation can also occur in the zygote cell, i.e. the initial cell formed when two gametes join.

A sporadic mutation can be the cause of an X-linked recessive disease (whereas it is unlikely for an autosomal recessive disorder) as a single mutation is enough in males to cause the disease. Males can be born affected due to a spontaneous gene mutation as a single abnormal gene copy is enough for the disease to become manifest; females can also be born carriers owing to random mutations.

Gametes are haployd cells, i.e., they only have one copy of the sex chromosomes (X or Y). To model the effects of spontaneous mutation on the diseases' distribution among the population, it is necessary to quantify the probability of inheriting an abnormal gene copy due to spontaneous mutation. Mutation may occur in any gene of the X chromosome, not all of them will cause a disease. From here on we consider only those mutations that occur on disease related genes. If the effects of de novo mutations are not considered, a carrier female has half of her gametes carrying a defective gene copy (X$^a$) and half hosting a normal gene (X$^A$) while the gametes of a healthy female contains only normal gene copies. All X gametes of an affected male carry the abnormal gene copy.

We let $\gamma$ denote the probability of a relevant mutation in the X chromosone of a gamete changing a normal gene copy into an abnormal one. Owing to spontaneous mutation, a healthy male has a probability equal to $\gamma$ of transmitting a defective gene copy to his daughters; consequently, $1-\gamma$ is the probability of transmitting a normal gene copy. In contrast, all the gametes of an affected male host only defective X chromosomes; hence the daughter of an affected male will receive a defective X chromosome with probability one.

Similarly, a healthy woman has a probability equal to $\gamma$ or $1-\gamma$ to transmit a defective gene copy or a normal gene copy to any of her progeny respectively. The combinations of these probabilities merged with the female and male X-linked inheritance pattern lead to the probabilities of having a healthy, carrier or affected child reported in Tables 1 and 2. Specialized medical literature, (see for example [12,37]) report $\gamma$ ranges in $[0, 10^{-5}]$; the rate of de novo mutations varies widely among different geographic regions, and depends on a number of factors, including environmental exposure to mutagenic agents, breadth of gene sequence and ability of the cell machinery to actually repair or correct the mutations [13].

Recent studies have shown that the longer is the disease's responsible gene the greater is the probability of spontaneous mutations.

Finally we notice that some recent studies suggest different mutation rates between females and males. This topic is currently discussed in the scientific community and there are not shared opinions (see for instance [42] [39]); the above differentiation is beyond the purpose of our work which specifically refers to diseases' responsible gene; hence in what follows we will adopt the mutation rate $\gamma$ to both females and males changing the non defective gene $A$ into $a$.

---

2.2. Individual fitness

The contribution of an individual to the next generation varies according to the individual's health and it is usually referred as individual's fitness [30,33]. Affected individuals may be in various ways disadvantaged to reproduce compared with healthy ones: this might be due to a lower likelihood to survive, hence to reach reproduction age, or to reduced fertility. These individuals, accordingly, will have relatively less children compared with healthy ones thus reducing the probability of passing the disease responsible gene to the next generation. This selective difference mechanism leads to a gradual evolution in the population's proportions among healthy, affected and carriers.

Following [32], the fitness factors $w_i$ are the product of the viability, i.e., the probability that an individual survives until the reproduction age, and the fertility, i.e., the average number of individual progeny. Hence the fitness factor $w_{i}$-also called fertility factor-models reproduction capacity of individuals of class $i$ according to their health conditions. The parameters $w_{i}$ are positive and are bounded by clinical considerations: the more severe the disease the smaller are their values in affected individuals.

---

3. Mathematical model of X-linked recessive diseases

In an attempt to further improve the predictive value of genetic disease' spread within a population, we designed a detailed model of the epidemiology of X-linked recessive diseases. Our model fits in the category of structured models [16,17]. The studied population is divided into homogeneous groups in relation to some major parameters (such as subject's age, sex or health conditions).

We divide the population into five classes, according to sex and health conditions (namely healthy and affected men, healthy, carrier and affected women) and, unlike other work in the literature, we choose the number of individuals in each class as system variables. Thus we introduce variables $x_1$ and $x_2$ to denote the number of healthy and affected males respectively and $x_3$, $x_4$ and $x_5$ to denote the number of healthy, carrier and affected females. This allows to reproduce the population distribution in the five classes at each generation and to predict the number of healthy, affected and carrier individuals at time.

We make the following assumptions:

• time is discrete and denotes generations;

• generations are non-overlapping thus each person breeds with a person of the same generation;

• the population is constant and equal to $2N$, that is,

$x_1(k)+x_2(k)+x_3(k)+x_4(k)+x_5(k)=2N \; \forall k;$

(1)

• in each generation there is an equal number of males and females, that is,

$x_1(k)+x_2(k)=x_3(k)+x_4(k)+x_5(k) \; \forall k.$

(2)

• The number of sons (equal to number of daughters) of each couple varies according to parents' health conditions and this is modeled through the fitness factors $w_{i}$ for $i=1 \dots 5$ that will be illustrated below.

• Spontaneous genetic mutations are modeled according to their probability of occurrence described in Section 2.1.

The number of males born from a person of class $i \in \left\{1, 2\right\}$ with a person of class $j\in\left\{3, 4, 5\right\}$ is:

$\frac{1}{2}w_{i}w_j\frac{x_i}{(x_1+x_2)}\frac{x_j}{(x_3+x_4+x_5)}(x_1+x_2+x_3+x_4+x_5) \\ = \frac{w_iw_jx_ix_j}{N}$

(3)

where $w_{i}$ is the fitness factor of a person of class $i$.

For example consider a couple formed by a healthy father (an individual of class 1) and a healthy mother (an individual of class 3); $\gamma$ sons of this couple (on average) will be affected and $\tilde{\gamma}$ will be healthy. Consider all couples of class $1, 3$ in the population; the number of their healthy sons at the next generation will increase of

$\frac{\tilde{\gamma}}{2}w_{13}\frac{x_1}{(x_1+x_2)}\frac{x_3}{(x_3+x_4+x_5)}(x_1+x_2+x_3+x_4+x_5).$

(4)

In view of the assumptions (1) and (2) system dynamics can be written using one state variable of the male population (i.e., $x_1$ or $x_2$) and two variables of the female population (i.e., two among $x_3$, $x_4$ or $x_5$). In the following we will restrict our analysis to healthy males and healthy and carrier females dynamics.

Let $x=[x_1 \; x_3\; x_4]^T$ be the state variables vector. The population dynamics are described through the following nonlinear, discrete-time system:

$x^+=f(x)$

(5)

where $x^+(k):= x(k+1)$, $f=[f_1 \; f_3 \; f_4]^T$ and

$f_1(x) =\frac{1}{N} \Big[w_{13}\tilde{\gamma}x_1x_3+ \frac{1}{2}w_{14}\tilde{\gamma}x_1x_4 +w_{23}\tilde{\gamma} (N-x_1)x_3 \\ + \frac{1}{2}w_{24}\tilde{\gamma} (N-x_1)x_4\Big]$

(6a)

$f_3(x) =\frac{1}{N}\Big[w_{13}\tilde{\gamma}^2 x_1x_3+\frac{1}{2} w_{14}\tilde{\gamma}^2 x_1x_4\Big]$

(6b)

$f_4(x) =\frac{1}{N}\Big[2w_{13}\gamma\tilde{\gamma} x_1x_3+ w_{23}\tilde{\gamma}(N-x_1)x_3 \\ + w_{15}\tilde{\gamma}x_1(N-x_3-x_4) + \frac{1}{2}w_{24}\tilde{\gamma} (N-x_1)x_4 \\ +\frac{1}{2} w_{14}(1+\gamma-2\gamma^2) x_1x_4\Big]$

(6c)

with

$\tilde{\gamma} := 1 - \gamma \; \mbox{and} \; w_{ij}:=w_iw_j \, .$

Term (4) is the first element in the right side of equation (6a); all other terms are derived with similar reasoning.

System (5)-(6) has been designed according to the inheritance pattern of X-linked recessive disease described in Section 2, the spontaneous mutation occurrence in Section 2.1, and given the newborn children as in (3).

We explicitly note that negative solutions of the dynamic equations (5)-(6) are not meaningful for an epidemiological model. The positivity of the state trajectories of system (5)-(6) is guaranteed assuming $w_{ij}$ within appropriate sets that may be derived exploring systems equation; an example of how these sets may be inferred is reported in the Appendix 8.

For the analysis in the following sections, it is helpful to rewrite system (5)-(6) in the form reported below:

$x^+=Ax+x_1Bx$

(7)

with

$A=\begin{bmatrix} 0&a_{12} &a_{13} \\0 &0 &0 \\ a_{31}& a_{32}& a_{33} \end{bmatrix}, B=\begin{bmatrix} 0&b_{12} &b_{13} \\0 &b_{22} &b_{23} \\ 0& b_{32}& b_{33} \end{bmatrix}$

(8)

where:

$a_{12} =a_{32}=w_{23}\tilde{\gamma} \\ a_{13} =a_{33}=\frac{1}{2}w_{24}\tilde{\gamma} \\ a_{31} =w_{15}\tilde{\gamma} \\ b_{12} =\frac{1}{N}\tilde{\gamma}(w_{13}-w_{23})\\ b_{13} =\frac{1}{2N}\tilde{\gamma}(w_{14}-w_{24}) \\ b_{22} =\frac{1}{N}w_{13}\tilde{\gamma}^2 \\ b_{23} =\frac{1}{2N}w_{14}\tilde{\gamma}^2 \\ b_{32} =\frac{1}{N}\left[2w_{13}\gamma\tilde{\gamma}-w_{23}\tilde{\gamma}-w_{15}\tilde{\gamma} \right] \\ b_{33} =\frac{1}{N}\left[\frac{1}{2}w_{14}(1+\gamma-2\gamma^2)-w_{15}\tilde{\gamma}-\frac{1}{2}w_{24}\tilde{\gamma}\right].$

---

4. System's equilibrium analysis

System (7) may admit many equilibrium points; according to [38] less than $2^3$. One is the zero vector $\bar{x}=(0, 0, 0)$ which implies $x_2=x_5=N$; when a population reaches this equilibrium, individuals are all affected.

Other equilibrium states depend on $w_i, \; i=1, \dots, 5$, and $\gamma$; if one component of the equilibrium point is either negative or greater than $N$, then it is unfeasible due to model's hypothesis.

Let $\bar{x}$ be a generic equilibrium point for system (7) and define $y \triangleq x -\bar{x}$. We rewrite (7) as

$y^+ +\bar{x}=A(y+\bar{x})+(y_1+\bar{x}_1)B(y+\bar{x}).$

Since $\bar{x}=A\bar{x}+\bar{x}_1B\bar{x}$ we obtain

$y^+ =\tilde{A}y+y_1 B y$

(10)

where

$\tilde{A}=A+\begin{bmatrix} \bar{x}_3b_{12} +\bar{x}_4b_{13}&\bar{x}_1b_{12} &\bar{x}_1b_{13} \\ \bar{x}_3b_{22} +\bar{x}_4b_{23} &\bar{x}_1b_{22} &\bar{x}_1b_{23} \\ \bar{x}_3b_{32}+ \bar{x}_4b_{33}& \bar{x}_1b_{32}& \bar{x}_1b_{33} \end{bmatrix}.$

Note that convergence of system (10) to $\bar{y}=(0, 0, 0)$ corresponds to convergence of system (7) to $\bar{x}$.

Theorem 4.1. If matrix $\tilde{A}$ is Schur then, the equilibrium state $\bar{x}$ of system (10) is asymptotically stable with

$\Lambda = \{ x\, |\, (x-\bar{x})^TH(x-\bar{x})< r \}$

(11)

as a region of attraction where

$\begin{array}{ll} r = \delta^2/\lambda_4 & \delta : = \sqrt{\left(\frac{b}{\lambda_1}\right)^2 + \frac{\lambda_1}{\lambda_3}} -\frac{b}{\lambda_1} \\ \lambda_1 := \lambda_{\min}(H^{-1}G) & \lambda_2:= 1/\lambda_{\min}(H) \\ \lambda_3 := \lambda_{\max}(H^{-1}B^THB) & \lambda_4 := \lambda_{\max}(H^{-1}CC^T) \\ b:= \lambda_2 \| B^T H \tilde{A}\| & C=[1\, \, 0 \, \, 0]^T \end{array}$

and $G$ and $H$ are positive definite symmetric matrices satisfying

$\tilde{A}^TH\tilde{A} - H = - G.$

(12)

The proof is carried out considering the following Lyapunov function

$V(y)=y^T H y \; \mbox{with} \; y = x-\bar{x}.$

(13)

Detailed demonstration is reported in the Appendix 8.

Remark 1. Theorem 4.1 for equilibrium point $\bar{x}=(0, 0, 0)$ requires that matrix ${A}$ is Schur; ranges of the fertility factors that guarantee this condition are given in the Appendix 8.

---

4.1. Stability analysis of the origin

Convergence of system (7) to the equilibrium point $\bar{x}=(0, 0, 0)$ coincides with the epidemiological situation where all the individuals are affected. Since this is not desirable it is useful to know the domain of attraction of $\bar{x}=(0, 0, 0)$ and avoid it. As matter of fact any point in this region represents a population distribution among healthy, affected and carriers individuals that leads, towards generations, to all affected people. In what follows our gaol is to delineate a wider domain of attraction for $\bar{x}=(0, 0, 0)$ compared to $\Lambda$ in Theorem 4.1 and explore global stability. To this aim, we propose a candidate Lyapunov function for the zero equilibrium point of system (7); we will show through simulation that this Lyapunov function leads to a larger domain of attraction compared to $\Lambda$.

In addition, we give conditions for global asymptotic stability of $\bar{x}=(0, 0, 0)$ in terms of fertility factors $w_{ij}$.

Consider the following candidate Lyapunov function:

$V(x)=\alpha x_1+\beta x_3+\mu x_4$

(14)

with $\alpha$, $\beta$ and $\mu>0$. It holds:

$\Delta V =V(x^+)-V(x)=(\mu a_{31}-\alpha)x_1+\Big[(\alpha+\mu) a_{12} \\ -\beta\Big]x_3 +\Big[(\alpha+\mu) a_{13}-\mu\Big]x_4+ x_1 \Big[(\alpha b_{12}+\beta b_{22}\\ +\mu b_{32})x_3+(\alpha b_{13}+\beta b_{23}+\mu b_{33})x_4\Big].$

(15)

---

4.1.1. Local stability analysis

The sufficient conditions for local stability (i.e., on linear terms in (15)) are:

$\mu a_{31}-\alpha <0 \rightarrow w_{15} <\frac{\alpha}{\mu \tilde{\gamma}}$

(16a)

$(\alpha+\mu) a_{12}-\beta <0 \rightarrow w_{23} <\frac{\beta}{\tilde{\gamma}(\alpha + \mu)}$

(16b)

$(\alpha+\mu) a_{13}-\mu <0 \rightarrow w_{24} <\frac{2\mu}{\tilde{\gamma}(\alpha +\mu)}.$

(16c)

Let's assume that (16) hold; equation (15) can be rewritten as:

$\Delta V =(\mu a_{31}-\alpha)x_1 \\ +\Big[(\alpha+\mu) a_{12}-\beta + (\alpha b_{12}+\beta b_{22}+\mu b_{32}) x_1 \Big] x_3 \\ +\Big[(\alpha+\mu) a_{13}-\mu + (\alpha b_{13}+\beta b_{23}+\mu b_{33})x_1 \Big] x_4$

in order to have $\Delta V <0$, $x_1$ should verify the following conditions:

$\left\{ \begin{array}{l} x_1 <\displaystyle \frac{\beta- (\alpha+\mu) a_{12}}{ \alpha b_{12}+\beta b_{22}+\mu b_{32}}=\displaystyle \frac{k_1}{k_2} \\ x_1 <\displaystyle \frac{\mu-(\alpha+\mu) a_{13}}{\alpha b_{13}+\beta b_{23}+\mu b_{33}}=\displaystyle \frac{z_1}{z_2} \end{array} \right.$

with $k_1>0$ and $z_1>0$. Depending on values of $k_2$ and $z_2$ we can derive the domain of attraction $\Phi$ of the equilibrium point $(0, 0, 0)$:

(ⅰ) $k_2>0$, $z_2>0$:

$\Phi=\Big\{x\Big|x_1 <\text{min}\left(\frac{k_1}{k_2}, \frac{z_1}{z_2}, N\right)\Big\};$

(19)

(ⅱ) $k_2>0$, $z_2 <0$: since $z_1>0$ the domain of attraction is given by

$\Phi=\Big\{x\Big|x_1 <\text{min}\left(\frac{k_1}{k_2}, N\right)\Big\};$

(20)

(ⅲ) $k_2 <0$, $z_2>0$: since $k_1>0$ the domain of attraction is identified by

$\Phi=\Big\{x\Big|x_1 <\text{min}\left(\frac{z_1}{z_2}, N\right)\Big\};$

(21)

(ⅳ) $k_2 <0$, $z_2 <0$: the equilibrium point is globally asymptotically stable.

Conditions on $k_2$ and $z_2$ in $(i)$, $(ii)$, $(iii)$ or $(iv)$ depends on $w_{13}$ and $w_{14}$ values; for example, holding (16), $k_2>0$ and $z_2>0$ if both of the following two conditions apply:

$w_{13} >\frac{\alpha + \beta}{\tilde{\gamma}(\alpha+\beta \tilde{\gamma}+2\mu \gamma)}\\ w_{14} >\frac{2(\alpha + \mu)}{\alpha\tilde{\gamma}+\beta \tilde{\gamma}^2+\mu(1+\gamma-2\gamma^2)}.$

Proposition 1. Assume that $w_{ij}$ verify (16) and one of conditions in $(i)$, $(ii)$ or $(iii)$ holds, then the equilibrium point $\bar{x}=(0, 0, 0)$ of system (7) is asymptotically stable within $\Phi$.

---

4.1.2. Global stability analysis

Proposition 2. If $w_{ij}$ are small enough such that constraints (16) hold and

$w_{13} <\frac{(\alpha+\mu)w_{23}+\mu w_{15}}{\alpha + \beta \tilde{\gamma}+2\mu \gamma}$

(23a)

$w_{14} <\frac{[(\alpha+\mu)w_{24}+2\mu w_{15}]\tilde{\gamma}}{\alpha \tilde{\gamma}+ \beta \tilde{\gamma}^2+\mu(1+ \gamma-2\gamma^2)},$

(23b)

then system (7) is globally asymptotically stable about $(0, 0, 0)$.

Proof. Constraints (23) on $w_{ij}$ assure that $k_2 <0$, $z_2 <0$.

Remark 2. So far we have introduced two Lyapunov functions; the first in (13) is suitable to analyze convergency properties of any equilibrium in system (10)-that is equivalent to system (5)-(6)-and gives a domain of attraction for the analyzed point. After that we designed the Lyapunov function (14) to specifically analyze the zero equilibrium point convergency properties. Compared to (13) the Lyapunov function (14) allows to:

1. Obtain the sufficient conditions for system's global asymptotic stability referred to $\bar{x}=(0, 0, 0)$.

2. Avoid conservative relations, such as norm and eigenvalue inequalities, for determining the domain of attraction; this could imply that region $\Phi$ is larger than $\Lambda$ although we were not able to formally prove this.

3. Relax the conditions of local stability for two state variables, namely $x_3$ and $x_4$.

---

4.2. Choosing optimal values of the parameters $\alpha$, $\beta$ and $\mu$

In this section we present a procedure to select the parameters' value $\alpha, \beta, \mu$ of the Lyapunov function in (14). Due to the critical role of the zero equilibrium point, we choose these values such that the domain of attraction $\Phi$ of the equilibrium point (0, 0, 0) is maximized.

The optimal values $\alpha^*$, $\beta^*$ and $\mu^*$ of the parameters are obtained by solving the following optimization problem:

$\underset{\alpha, \beta, \mu}{\operatorname{max}} \, \Phi(x, \alpha, \beta, \mu)$

subject to

$\mu a_{31}-\alpha <0 \\ (\alpha+\mu) a_{12}-\beta <0 \\ (\alpha+\mu) a_{13}-\mu <0 \\ \alpha, \beta, \mu >0.$

The objective function to be maximized is one of the regions of attraction (19), (20) and (21) according to $k_2$ and $z_2$ signs; the constraints are given by the sufficient conditions (16) for the local stability of the zero equilibrium point.

Genetic Algorithm (GA) method has been exploited for solving this constrained nonlinear optimization problem which is not well suited for standard optimization methods. GA peculiar mechanisms to avoid local minima make the use of a GA reasonable in facing problems with complex nonlinear objective function (like in our case) in place of gradient based methods which are designed to be fast and efficient for finding local minima.

---

5. Simulations

In what follows we present some simulations of system (7) that is equivalent to (5)-(6). We show results obtained in Section 4 on the local stability of the equilibrium point $\bar{x}=(0, 0, 0)$; namely we evaluate the domains of attraction $\Lambda$ and $\Phi$ obtained through the Lyapunov functions (13) and (14) respectively. Moreover we give results related to a non zero equilibrium and show its convergency according to Theorem 4.1. Finally we depict the trajectories for some of the analyzed scenarios.

---

5.1. Domain of attraction for $\bar{x}=(0, 0, 0)$.

We assign three different sets of values to $w_i$ parameters, reported in Table 3; they have been chosen such that matrix $A$ in system (7) is Schur (i.e., conditions (27) are satisfied).

Using Theorem 4.1 we determine the domain of attraction $\Lambda$ of the equilibrium point $\bar{x}=(0, 0, 0)$ according to (11); $\Lambda$ is an ellipsoid centered in the equilibrium point with semi-principal axes equal to the square roots of the reciprocals of matrix $H$ eigenvalues multiplied by $r$. Matrix $H$ is the solution of discrete Lyapunov equation (12) where matrix $G$ has been chosen equal to the identity matrix. Figure 2 depicts the domains of attraction $\Lambda$ for each scenario in Table 3, where the corresponding values of $r$ are reported.

To apply results in Proposition 1 we evaluate the coefficients of Lyapunov function (14) for each scenario in Table 3. To this aim we solve the optimization problem in Section 4.2 and we compute the domains of attraction $\Phi=\{x|x_1<\underline{x}_1\}$; the obtained $\alpha^*, \beta^*, \mu^*$ and $\underline{x}_1$ for the three scenarios are reported in Table 4.

Figure 3 depicts the domains of attraction $\Phi$ compared to the ellipsoids $\Lambda$. Note that in each of the reported scenarios $\Phi$ is sensibly greater than $\Lambda$; indeed for the first two scenarios $\Phi$ corresponds to system states space.

---

5.2. Equilibrium point other than $\bar{x}=(0, 0, 0)$.

We consider $w_1=0.9$, $w_2=2$, $w_3=1$, $w_4=1$ and $w_5=0.5$, $\gamma=10^{-4}$ and $N=150$; for this choice of parameters system (5)-(6) has the following two equilibrium points: $\bar{x}_A=(0, 0, 0)$, $\bar{x}_B=(72.5, 26.5, 68.9)$. Note that these parameters do not satisfy (27); thus matrix A in system (7) is unstable, and system (5)-(6) will not converge to $\bar{x}_A$. Finally matrix $\tilde{A}$ in system (10) evaluated at $\bar{x}_B$ is stable, and Theorem 4.1 applies. The corresponding domain of attraction $\Lambda$ has $r=341.45$ and it is depicted in Figure 4.

---

5.3. State portrait

Figure 5 depicts trajectories of system (5)-(6) in state space with the $w_i$ parameters assigned according to the scenario 1 in Table 3. The initial conditions are chosen within the corresponding domain of attraction shown in Figure 2; as expected, all trajectories asymptotically go to zero.

Figure 6 represents trajectories of system (5)-(6) and convergency to $\bar{x}_B$ when initial conditions are inside the domain of attraction obtained in Section 5.2.

---

6. Sensitivity analysis

Local sensitivities have been performed to analyze the dependency of the dimension of the domain of attraction $\Lambda$ in (11) on system's fitness factors. Local sensitivity analysis is performed in a static way; the output obtained changing each parameter in a predefined interval and with constant step are evaluated. Namely we quantify the influence of couples' fitness factors $w_{ij}$ on $r$ in (11). The analysis is performed perturbing one parameter at time on M different values equally distributed on an interval. We set $w_{ij}=\bar{w}_{ij}$, and obtain the nominal trajectory $r^{\rm{nom}}$; $r^{\rm{pert}}$ is obtained perturbing one parameter from its nominal value. For each perturbed parameter, we compute the following function:

$\Delta E_l= \left( r_{l}^{\rm{pert}}-r^{\rm{nom}}\right)^2,$

where $l=1, \ldots, M$ denotes the values of the considered perturbed parameter.

To compute the nominal trajectory we consider $\gamma=10^{-4}$, $N=150$ and we fix the nominal conditions $[\bar{w}_{13}, \bar{w}_{14}, \bar{w}_{15}, \bar{w}_{23}, \bar{w}_{24}]=[0.5, 0.5, 0.5, 0.5, 0.5]$. The perturbed trajectories are obtained letting $w_{ij}$ vary in the range $[0, 1.3]$, with a constant step of $0.01$. Only one parameter varies per time; all others are set at their nominal value. To each combination of $w_{ij}$ values corresponds a set of system's equilibrium points that always include $\bar{x}=(0, 0, 0)$; in order to perform a meaningful comparison between the domains of attraction we only referred to the $\Lambda$ of the zero equilibrium. The range in which the $w_{ij}$ vary assures that the matrix $A$ is always Schur.

Figure 7 depicts the sensitivity functions with respect to each $w_{ij}$. A flat sensitivity curve represents a little influence of the parameter on the output; a steep curve implies a great influence. Function $r$ in (11) appears to be most sensitive to $w_{15}$ and then to $w_{24}$, that are related to the sufficient conditions for the local stability of the system (see equation (27)).

---

7. Discussion

Genetic epidemiological studies have the ultimate goal of estimating the spread of a disease in a population. This is quite a recent and rapidly expanding research filed, but the implications on individual or population health are still unclear [36]. Findings in genetic epidemiology are exploited to design genetic screening campaign aiming at assessing, through DNA sequencing tests, how many people are affected or carry the disease's responsible gene. Usually females carring a defective gene copy don't show any symptoms of the disease, thus their number is usually hard to estimate. State $x_4$ of model (5)-(6) predicts the number of carrier in a population and may be successfully applied to genetic epidemiology.

Genetic epidemiology may also contribute to measure the effect of environment's risk factors on disease spread. Urban areas exposed to toxins, radiations or other agents are periodically examined; sampled data on the spread of genetic diseases are compared with predicted values determined through mathematical models. By simulating model (5)-(6) one can evaluate the effects of fitness factors variations on spread of the disease. Simulations results may help in taking decisions on appropriate corrective strategies aiming at preventing the spread of the disease.

Finally technological and other advances will allow the potential of genetic epidemiology to be revealed over the next few years, and the establishment of large population-based resources for such studies (biobanks) should contribute to this endeavor.

Consider an X-linked disease in an assigned population $2N$; parameters in model (5)-(6)-the mutation rate $\gamma$ of the disease responsible gene and the fertility factors $w_i \; i=1, \dots, 5$-can be derived through clinical considerations or found in specialized medical literature. Applying the results on global asymptotic stability of equilibrium $\bar{x}=(0, 0, 0)$, one can conclude that the population will asymptotically converge to all individuals affected if the assigned values of $w_{i}$ verify (16) and are within the range defined in (23). For the other equilibrium points of system (5)-(6) one can verify if matrix $\tilde{A}$ is stable and derive the corresponding region of attraction using Theorem 4.1 results. In what follows we apply our model and the results on stability analysis to two X-linked recessive diseases: red and green color blindness and Hemophilia A.

---

7.1. Application 1

Color blindness, or color vision deficiency, is the inability or decreased ability to see color, or perceive color differences, under normal lighting conditions. Color blindness affects a significant percentage of the population. The most common cause of the disease is a fault in the development of one or more sets of retinal cones that perceive color in light and transmit that information to the optic nerve. This type of color blindness is usually an X-linked condition as genes that produce photopigments are carried on the X chromosome; if some of these genes are missing or damaged, color blindness will be expressed in males with a higher probability because in females a functional gene on only one of the two X chromosomes is sufficient to yield the needed photopigments. According to statistical datas, [23] and [1], color blindness affects a significant number of people, although exact proportions vary among groups. Isolated communities with a restricted gene pool sometimes produce high proportions of color blindness, including the less usual types. Examples include rural Finland, Hungary, and some of the Scottish Isles. We apply our model to study the red and green color blindness diffusion in Scottish Isles. As reported in the National Records of Scotland, in 2011 Scottish Isles had a population of 87252. According to data in [34], the mutation rate for red and green color blindness ($\gamma$ in our model), is approximately $2\cdot 10^{-5}$. Red and green color blindness does not effect individual's ability to survive or procreate thus unitary fitness factors ($w_i=1, \; i=1, \dots, 5$) can be assigned. Substituting the above parameters' values in system (5)-(6) the following equilibrium points are found: $\bar{x}_A=(0, 0, 0)$, $\bar{x}_B=(40231.5,37100.5,6263.7)$. Matrix A in system (7) is not Schur, thus model (5)-(6) will not converge to equilibrium $\bar{x}_A$, that is the population will never converge to the state where all people are affected. Moreover matrix $\tilde{A}$ in (10) is Schur if evaluated at $\bar{x}_B$. Applying Theorem (4.1), the domain of attraction $\Lambda$ for equilibrium $\bar{x}_B$ is $\Lambda = \{ x\, |\, (x-\bar{x})^TH(x-\bar{x}) < 1.23\cdot10^{6}\}$. Interestingly enough the distribution among healthy, affected and carrier people in equilibrium $\bar{x}_B$ is very close to the one in literature [1].

---

7.2. Application 2

Hemophilia A disease, is a hereditary bleeding disorder caused by a lack of blood clotting factor Ⅷ, a protein encoded by FⅧ gene placed on the X chromosome. It is largely an inherited disorder, that is the spontaneous gene mutation rate of the diseases ($\gamma$ in our model) can be considered negligible [6]; therefore we can set $\gamma=0$. Affected males show a reduced reproduction capacity related to the severity of the disease symptoms; carrier females do not usually show any sign of the disease [9]. Due to the severity of the disease, one can assign fertility factors equal to zero to affected individuals ($w_2=w_5=0$); based on clinical observations the fertilities factors of the other classes can be assigned as follows: $w_1=w_3=1$ and $w_4=0.6$. We studied the diffusion of the hemophilia A in the Italian population ($2N$ in our model) that was equal to 59685228 individuals in 2012 (see [20]). Substituting the above parameters' values in system (5)-(6), the following equilibrium points are found: $\bar{x}_A=(0, 0, 0)$, $\bar{x}_B=(N, N, 0)$. The domain of attraction of equilibrium $\bar{x}_A=(0, 0, 0)$ can be derived applying results in Proposition 1. We compute the coefficients of Lyapunov function (14) by solving the optimization problem in Section 4.2 through GA, obtaining the optimal values $\alpha^*=0.01$, $\beta^*=4.83$ and $\mu^*=2.51$ and we compute the domain of attraction $\Phi=\{x\Big|x_1<5.95*10^7\}$. Indeed for almost any initial population distribution system (5)-(6) will converge to $\bar{x}_A$; different fertility factors datas, gathered from direct population analysis, could lead to a different population distribution.

---

8. Conclusions

We developed a discrete time nonlinear model for X-linked recessive diseases aiming at describing the spread of such diseases in a population. The model accounts for de novo mutations on the inheritance pattern and distinct fitness factors. Under the assumption of constant population size we analyzed system's properties and performed stability analysis of equilibrium points through Lyapunov second method.

Extensions of the present work should consider sensitivity analysis to parameters variations as well as the effect on the results when weakening some of the model assumptions such as the hypothesis that the population has a constant size or an equal number of females and males.

Future developments will assign different de novo gene mutation rates to females and males as well as include de novo mutation possibly changing a defective gene into a normal one.

---

Appendix

Proof of Theorem 4.1.

Proof. Consider the candidate Lyapunov function:

$V(y)=y^T H y \mbox{with} y = x-\bar{x}.$

Then,

$\Delta V(y) =V(y^+)-V(y) \\ =[\tilde{A}y+y_1By]^T H [\tilde{A}y+y_1By]-y^T H y \\ =y^T (\tilde{A}^T H \tilde{A}-H)y +2 y^TB^T H \tilde{A}y^T y_1 + y^TB^T HByy_1^2 \\ = - y^T Gy + 2 y^TB^T H \tilde{A} y y_1 + y^TB^THBy y_1^2 \\ \le -y^T Gy +2\|y\|^2 \| B^T H \tilde{A}\| |y_1| + y^TB^THBy |y_1|^2.$

(25)

Since $H^T = H >0$, we can use the following as upper bounds for terms in inequality (25):

$\begin{array}{lll} -y^T Gy\leq - \lambda_1 y^THy & \lambda_1 := \lambda_{\min}(H^{-1}G) \\ \|y\|^2 \le \lambda_2 y^THy & \lambda_2 := 1/\lambda_{\min}(H) \\ y^TB^THBy \leq \lambda_3 y^THy & \lambda_3 := \lambda_{\max}(H^{-1}B^THB). \\ \end{array}$

Substituting these inequalities into (25) leads to:

$\label{eq:DeltaV2} \Delta V \le -V(y)( \lambda_1 -2b |y_1| - \lambda_3 |y_1|^2)$

(26)

where $b:= \lambda_2 \| B^T H \tilde{A}\|$. One may verify that $\lambda_1 -2b |y_1| - \lambda_3 |y_1|^2 >0$ if and only if $|y_1| < \delta$ where

$\delta : = \sqrt{\left(\frac{b}{\lambda_1}\right)^2 + \frac{\lambda_1}{\lambda_3}} -\frac{b}{\lambda_1} >0.$

Since $|y_1|^2 = (Cy)^TCy$ where $C=[1\, \, 0 \, \, 0]^T$ and $H^T = H >0$, exploiting properties on eigenvalues of symmetric and positive definite matrices we have

$|y_1|^2 \le \lambda_4 y^THy \lambda_4 := \lambda_{\max}(H^{-1}CC^T).$

Hence, whenever $V(y) < r = \delta^2/\lambda_4$ we have $|y_1| < \delta$. It now follows from (26) that

$\Delta V(y) < 0 \; \mbox{when} \; 0 < V(y) < r.$

Conditions for matrix A to be Schur. When $\bar{x}=\left(0, 0, 0 \right)$, matrix $\tilde{A}$ in system (10) is equal to matrix $A$ in (8).

Matrix $A$ in (8) is Schur iff :

$\tilde{\gamma} w_{24} + \tilde{\gamma}^2w_{15}w_{24} < 2.$

(27)

Proof. The characteristic polynomial of $A$ is given by

$p(\lambda)=\det(\lambda I -A) =\lambda^3-a_{13}\lambda^2-a_{13}a_{31}\lambda.$

Since $a_{13}, a_{31} \ge 0$, the criterion for a Schur polynomial reduces to

$a_{13}a_{31} < 1 \\ a_{13} < 1-a_{13}a_{31}$

which are equivalent to

$a_{13} +a_{13}a_{31} xy 1$

that is,

$\tilde{\gamma} w_{24} + \tilde{\gamma}^2w_{15}w_{24} < 2.$

Conditions on $w_{ij}$ for positiveness of system (5)-(6). In what follows we prove that system (5)-(6) is positive if the system is initialized in a positive state and if $w_{ij}$ are chosen within an appropriate set.

In particular we detail the proof for $f_1(x)$; analogous reasoning gives the conditions for the positiveness of $f_3(x)$ and $f_4(x)$.

Let's assume $w_{ij}\le 1$ $\forall i, j$, $w_1\ge w_2$, $w_3\ge w_4$, $\gamma\le 1$ and the initial states $0\le x_i(0)\le N$, $i=1, 3, 4$; if we rewrite (6a) as:

$f_1(x)=\frac{1}{N} \Big[w_{3}\tilde{\gamma}(w_1-w_2)x_1x_3 + \frac{1}{2}w_{4}\tilde{\gamma}(w_1-w_2) x_1x_4 +w_{23}\tilde{\gamma} Nx_3 +\frac{1}{2}w_{24}\tilde{\gamma}Nx_4\Big]$

we can deduce that, holding the previous assumptions, $x_1(1)\ge 0$. Consider now

$N-f_1(x) =N-\frac{\tilde{\gamma}}{2N}\Big[(2w_3 x_3 + w_4 x_4) \left(Nw_2+\left(w_1-w_2\right) x_1\right) \Big]\\ \ge N-\frac{\tilde{\gamma}}{2N}\Big[w_3(2 x_3 + x_4) \left(Nw_2+\left(w_1-w_2\right) x_1\right) \Big]\\ \ge N-\frac{\tilde{\gamma}}{2N}\Big[2 w_3 N \left(Nw_2+\left(w_1-w_2\right)N\right) \Big]\\ = N-\tilde{\gamma}w_{13}N$

(29)

Under the previous assumptions the condition $N-f_1(x)\ge 0$ holds, hence $x_1(1)\le N$. By mathematical induction, we prove $0\le x_1(k)\le N$, $\forall k>0$.

Note that conditions on $w_{ij}$ do not affect the comprehensiveness of our epidemiological model, since they simply require that affected individuals show a fitness factor less than or equal to healthy ones.

---

Acknowledgments

Authors are grateful to prof. Emanuele Durante Mangoni, MD PhD, from the Internal Medicine Unit, Second University of Naples, Italy, for helpful discussions and literature suggestions on X-linked diseases and their inheritance pattern and for his critical review of the manuscript.

---