International Journal of

1. Introduction

Infertility is referred to the inability of a woman to become pregnant after at least one year of regular sexual intercourse without using birth control. Infertility is seen in 10–15% of couples (1). The main causes of infertility include: ovulatory dysfunction, tuboperitoneal pathology, male factors, and uterine pathology (2). One of the most common causes of ovulatory dysfunction is the polycystic ovary syndrome. In 2003 Rotterdam, the presence of two out of three criteria for diagnosis of PCOS is essential. (i) Clinical and/or biochemical signs of hyperandrogenism, (ii) Oligomenorrhea or anovulation, and (iii) polycystic ovaries (with
the exclusion of related disorders) (3). In most cases, women with PCOS have reduced fertility due to low ovulation; therefore, the ovarian stimulation methods are used for a group of these patients due to their reduced ovarian response. Clomiphene Citrate is the major component of the ovarian stimulation treatment, causing 80–85% of the ovulation and leading to pregnancy in 40% of women (2) Gynecologists suggest that after a maximum of six months being in ovulation cycles with the use of Clomiphene Citrate, it is better for the patient to be treated with gonadotrophin. If the treatment cycles with Clomiphene Citrate and
Gonadotropin are determined and its therapeutic follow-up is controlled, and yet the pregnancy does not occur during the 9-12 treatment cycles, the ovulation stimulation method should not be relied upon. Therefore, using the Assisted Reproductive Techniques (ART), one of which is IVF, is recommended. Here, we present a model for women infertility with polycystic ovary syndrome, based on the individuals who went to the Yazd Infertility Research and Treatment Clinic. To this end, we show the total number of the infertile couples who at the time of t went to the clinic over the course of a year, with N(t). We have divided them into five groups, including: the newcomers who went to the clinic and their disease has not yet been confirmed, patients (those whose disease has been diagnosed to be and should be treated by this clinic), those treated with Clomiphene Citrate and Gonadotropin, those treated with ART, and improved patients
(those whose test was positive). We call this our model.
In the second section, we model infertility in women with PCOS using dynamic system methods and find the number of the secondary sufferers during the course of the disease that had returned to the disease, and we show it with R0.
In the third section, we obtain the equilibrium point of the model SIT1T2R, using geometric methods and examine its stability, asymptotically.
In the fourth section, we examine the stability of the disease asymptotically in the Q∗ equilibrium point.
In the fifth section, using the Stoke’s Theorem, we examine the non-intermittence of the therapeutic cycle.
In the sixth section, the model is solved numerically by using Rung-Kutta method, and the obtained data are summarized in the several graphs, which are analyzed in detail.

2. Materials and Methods

According to the previous section, we show the total number of the infertile women at time t with N(t). Therefore, N(t) = S(t) + I(t) + T1(t) + T2(t) + R(t).
Suppose 𝛼 and 𝛾 refer to the birth rate of the disease (the rate of the patient’s arrival at the clinic to diagnose and treat their disease) and the death rate (leaving the clinic without any result or not being treated), respectively. By ‘not being treated’, we mean those who were checked at the clinic, and the type of their treatment was diagnosed, but they discontinued their treatment. The rate of abortion and returning to the cycle of treatment is 𝛽. The treatment speed of those who have become pregnant through medication (clomiphene Citrate and gonadotropin) is 𝜂, and 𝜁 is the treatment speed of those who underwent IVF treatment, and their BHCG tests were positive.
We show the mean of the infertile people who had abortions or were with OHSS disease and returned to the treatment cycle at time t by 𝛽N(t). The probability of a disease returning in the susceptible woman is equal to  (4). Therefore, the number of new women patient per time unit is equal to 𝛽N  , and the number of infertile women who have the possibility of returning to the disease at time t is equal to 𝛽(𝜂T1 + 𝜁T2)S. D indicates the number of people who left the clinic without obtaining any result from the treatment or without being treated. The number of people prone to the disease entering group D is 𝛼S (Figure 1).



Figure 1. Representation of the susceptible couples.


Therefore, S′ = 𝛾N − 𝛼S − 𝛽(𝜂T1 + 𝜁T2)S.
We show the treatment rate in the patient group with 𝜆 and the treatment rate of 𝜆I who were under medical treatment with b. Therefore, the number of patients who used IVF treatment is equal to (1 - b)𝜆I. The number of people who left the clinic without any treatment is 𝛼I (Figure 2).



Figure 2. Patient couples are divided to three groups.


Hence I′ = 𝛽(𝜂T1 + 𝜁T2)S − (𝛼 + 𝜆)I.
We show the recovery rate in group T1 with 𝜇1. We divide 𝜇1T1 into two groups: the first group is the number of people who have recovered, and we show it with u1𝜇1T1, in which u1 is the recovery rate of 𝜇1T1 at time t. Therefore, (1 − u1)𝜇1T1 of these people enters T2 group. The number of people under medical treatment who left the clinic without any treatment is 𝛼T1 (Figure 3).



Figure 3. Infertile women under drug treatment.


Thus T′ 1 = 𝜆bI − (𝛼 + 𝜇1)T1.
We show the recovery rate in group T2 with 𝜇2. The recovery rate of 𝜇2T2 at time t is u2, therefore, the number of people treated with IVF is equal to u2𝜇2T2 and the number of people under treatment entering group D is equal to (1 − u2)𝜇2T2. The number of IVF patients who left the clinic without any result is 𝛼T2 (Figure 4).



Figure 4. Infertile women under IVF treatment.

According to Figure 4, we have: T′ 2 = 𝜆(1 − b)I + (1 − u1)𝜇1T1 − (𝛼 + 𝜇2)T2.
The number of susceptible people who left the group due to death is 𝛼R (Figure 5).



Figure 5. Cured infertile couples.

Hence R′ = u1𝜇1T1 + u2𝜇2T2 − 𝛼R.
The number of infertile women who used cure methods is denoted by 𝛾N, and those entered susceptible at t time by (1 − u2)𝜇2T2 and 𝛼N denotes the number of infertile women who leave the population at t time (Figure 6).



Figure 6. All of the infertile couples.


According to the aforementioned cases, the dynamic model of infertile patients with PCOS can be modeled as follows:

S′ = 𝛾N − 𝛼S − 𝛽S (𝜂T1 + 𝜁T2)
I′ = −𝜆I − 𝛼I + 𝛽S (𝜂T1 + 𝜁T2)
T₁′ = −𝜇1T1 − 𝛼T1 + 𝜆Ib
T2′= −𝜇2T2 − 𝛼T2 + (1 − u1) 𝜇1T1 + (1 − b)𝜆I
R′ = −𝛼R + u1𝜇1T1 + u2𝜇2T2
N′ = 𝛾N − (1 − u2) 𝜇2T2 − 𝛼N.

We assume that the rate of return to the disease in each infertile women is the constant value of 𝛽, and the number of secondary sufferers is equal to the number of infertile women who were treated, but returned to this disease afterward during the course of treatment. At first, we assume that the number of people in the population is equal to the number of people susceptible to the disease. To find R0, we replace N(t), and S(t) with S0 and deduce the following
system.



This linear system is divided into two parts. The first matrix is shown by F and is called returning matrix and the second matrix by K, which is called affected matrix. Thus, if






                                                                                               .


2.1. Analysis of women’s infertility based on asymptotical stability

Based on the biological hypothesis, we can deduce SIT1T2R model.
Prior to the disease, the free equilibrium point is Q0 = (S0, 0, 0, 0, S0), and we have the following result (6).

Theorem 2.1.1. If R0 < 1 then:
a) If 𝛾 −𝛼 > 0, then Q0 is the free equilibrium point which is locally stable but it is not asymptotically stable.
b) If 𝛾 − 𝛼 = 0, then the free equilibrium point is unstable.
c) If 𝛾 − 𝛼 < 0, then Q0 is locally asymptotically stable.
Proof. (a) The linearized matrix of model SIT1T2R at the point Q0 is


                                                                                     

The characteristic equation of this matrix is:




The roots of this equation are:                                                          







Then, for other roots equation (3.3) using Routh Hurwitz methods (7). As all of the parameters are positive, a1 = (3𝛼 + 𝜇2) + 𝜆 + 𝜇1 > 0. Moreover, a0 = 1 > 0.
We take a2 = [−b𝜆𝛽𝜂S0 + (𝛼 + 𝜆)(𝛼 + 𝜇1) + (𝛼 +𝜇2)(𝜆 + 𝛼) + (𝛼 + 𝜇1)(𝛼 + 𝜇2) − 𝛽𝜁S0(1 − b)𝜆].                                          
Additionally, a3 = [−b𝜆𝛽𝜂S0(𝛼 + 𝜇2) − b𝜆𝛽S0𝜁𝜇1(1 − u1) + (𝛼 + 𝜇1)(𝛼 + 𝜇2)(𝛼 + 𝜆) −  (𝛼 + 𝜇1)𝛽𝜁S0(1 − b)𝜆].
Since R0 < 1, then
           

                                                                                                        

Once again as R0 < 1, then (𝛼+𝜆)(𝛼+𝜇1)(𝛼+𝜇2) > 𝛽𝜁S0(1−u1)𝜇1b𝜆+𝛽𝜂S0b𝜆(𝛼+𝜇2)+𝛽𝜁S0(1−b)𝜆2(𝛼+𝜇1).

Hence,

                                                                                                           



Since







Therefore, a1a2 − a0a3 = (𝛼 + 𝜇1)2(𝜆 + 𝛼) + (𝛼 +𝜆)2(𝜇2 + 𝛼) + (𝛼 + 𝜇2)2(𝛼 + 𝜆) + (𝛼 + 𝜆)2(𝜇1 + 𝛼) + (𝛼 +𝜇2)(𝛼 + 𝜆)(𝛼 + 𝜇1) + (𝛼 + 𝜇1)2(𝛼 + 𝜇2) + 𝑏𝜆𝛽𝜂s0(𝛼 +𝛼) − b𝜆𝛽𝜂S0(𝛼 + 𝜇1) − 𝛽S0𝜁𝜆2(1 − b)(𝜇1 + 𝛼) − (𝛼 +𝛼)𝛽𝜁S0𝜆(1−b)+𝛽S0𝜁𝜆b(1−u1)𝜇1 +𝛽S0𝜂𝜆b(𝜇2 +𝛼) >0 ⇒ Δ2 > 0.


                                                                                                         
By using (3.5) and (3.6), we obtain


                                                                                                             
Hence, by Hurwitz criterion, the real parts of the roots of (3.3) are all negative. If (3.2) is positive, then Q0 is not locally asymptotically stable.

(b) If (3.2) is zero, then Q0 is unstable.
(c) If 𝜆′₂ = 𝛾 − 𝛼 < 0, then the free equilibrium point Q0 is locally and asymptotically stable. (8, 9)
If 𝛾 > 𝛼, the number of patients who leave the clinic without result is less than the number of patients who enter the clinic. This demonstrates that considered treatment is useful. For a number of patients with Polycystic Ovary function, 𝛾 > 𝛼 means local stability of free point; that is not asymptotically stable. If 𝛾 < 𝛼, it shows the asymptotical stability of system, and it states that patients who come to the center have an acute infertility problem. If 𝛾 = 𝛼, system is unstable,
it means all the patients who come to the center leave there without any result. The evaulation of the treatment process depends on these three cases.







Then, Q∗ = (S∗, I∗, T1∗ , T2∗ , N∗) shows the other equilibrium point. Since the Q∗ point is located in the interior of the positive space, we must have:

                                                                         



2.2. Sign stability of the dynamicalmodel of women’s infertility

When we model a biological phenomenon with differential equations, the values of the parameters used in them are error-prone. If we deal with a square matrix, then it is highly important to know how much the stability of this matrix depends on the elements of that matrix and how sensitive is it to the variations of the matrix elements. Therefore, in this section, using the following definitions and Theorems, we examine the sign stability of the system (1) at the equilibrium point of Q∗. To state the Theorem (2.2.3.), we use some concepts of graph theory as follows.

Definition 2.2.1.(10): An n × n square matrix A = [aij] is said to be sign stable if every n×n square matrix B = [bij] of the same sign pattern (i.e., signbij = signaij for all i; j = 1, 2, ..., n) is a stable matrix.

For an n×n square matrix A = [aij], we can obtain an undirected graph GA, whose vertex set is V = 1,2,..., n and edges are {(i, j) ∶ i ≠ j; aij ≠ 0 ≠aji; i, j = 1, 2, ..., n}. Also, a directed graph DA can also attach to A with the same vertex set and edges {(i, j) ∶ i ≠ j; aij ≠ 0 ≠ aji; i, j = 1, 2, ..., n}. A kcycle of DA is a set of distinct edges of DA of the form: {(i1, i2), (i2, i3), ...., (ik−1, ik), (ik, i1)}: Let RA = {i ∶aii ≠ 0} ⊂ V, which are the numbers for them the corresponding element in the main diagonal of the matrix is not zero. An RA-coloring of GA is a partition of its vertices into two sets, black and white (one of which may be empty). Such that each vertex in RA is black, no black vertex has precisely one white neighbor, and each white vertex has at least one white neighbor. A V − RA complete matching is a set M of pairwise disjoint edges of GA such that the set of vertices of the edges in M contains every vertex in V − RA. By applying this concepts, we are now able to state the following Theorem.
Theorem 2.2.2.(10): An n × n real matrix A = [aij] is sign stable if it satisfies the following conditions:
(i) aij ≤ 0 for all i, j;
(ii) aijaji ≤ 0 for all i ≠ j;
(iii) The directed graph DA has no k-cycle for ≥ 3;
(iv) In every RA-coloring of the undirected graph GA all vertices are black; and
(v) The undirected graph GA admits a RA complete matching.

Theorem 2.2.3. If 𝛾 > 𝛼, then the matrix model SIT1T2R at the equilibrium Q∗ is not sign stable.

Proof. The matrix of the linearize system model of SIT1T2R at the equilibrium Q∗ is given by





Where P∗ = I∗. Theorem (2.2.2.) requires





Which is not sign stable. Since a55 > 0, then Theorem 2.2.2. implies that the matrix A is not sign stable.

Theorem 2.2.3. implies that in SIT1T2𝑅 model at equilibrium point Q∗ is not sign stable. This is a

Theorem 2.2.2.(10): An n × n real matrix A = [aij] is sign stable if it satisfies the following conditions:
(i) aij ≤ 0 for all i, j;
(ii) aijaji ≤ 0 for all i ≠ j;
(iii) The directed graph DA has no k-cycle for k ≥ 3;
(iv) In every RA-coloring of the undirected graph GA all vertices are black; and
(v) The undirected graph GA admits a RA complete matching.

Theorem 2.2.3. If 𝛾 > 𝛼, then the matrix model SIT1T2R at the equilibrium Q∗ is not sign stable.

Proof. The matrix of the linearize system model of SIT1T2R at the equilibrium Q∗ is given by positive case from the perspective of the disease treatment, since the sign stability of SIT1T2R model is equivalent to the return of disease, which has not been reversed for a long time; as a result, the problem of infertility and reproduction still remains for a high percentage of the population.

2.3. Reversibility of the infertility cycle with the PCOS factor in an SIT1T2R model is not intermittent

The infertility treatment with PCOS initiates as a course of treatment with Clomiphene Citrate and Gonadotropin in some menstrual cycles, which, by the failure of medical treatment, changes the treatment with the IVF therapy cycle. These patients will stay in that therapy course until they achieve the desirable result from it. From a medical point of view, this therapeutic cycle will not include the return to medical treatment during its treatment course. In this section, we
prove that in system (1), the return of infertility is not circular in the cycle of the disease with the PCOS factor. In other words, the disease has no circulation.
First, we state the Stoke’s Theorem.

Theorem 2.3.1. (Stoke’s Theorem) (9): Suppose that M is a bordered, oriented, compact, k dimensional Manifold, and also suppose that 𝜔 is a smooth k − 1 form on M. In this case:




Where 𝜕M, as described earlier, is an oriented border.

Now, with the help of Stoke’s Theorem, we state the Theorem on the closed circuits in this model.

Theorem 2.3.2. The SIT1T2R model does not have a closed circuit.

Proof. Assumption by contradiction: If there is a closed circuit like C with the rotation period of T, then we consider the Manifold M as a region in R6 with the border C , and assume 𝜔 as a 5-form on M, which is defined as follows:




Hence, ∫M d𝜔 < 0 and ∫𝜕M=C 𝜔 = 0. However, it is concluded from Stoke’s Theorem that 0 = ∫M d𝜔 < ∫𝜕M=C 𝜔 = 0 and this is a contradiction. Therefore, the system (1) has no closed circuit; in other words, the reversibility of the therapeutic cycle of this disease is not intermittent.

2.4. Numerical solutions of the model 𝑆𝐼𝑇 1𝑇 2𝑅

In this section, we solve SIT1T2R model numerically via Rung-Kutta method (11). First, we define the following functions:







Hence, S0, I0, (T1)0 , (T2)0, N0, and R0 are the initial conditions.

We take h = b−a M , where M > 0 and tj = a + jk for j = 0, 1, 2, … , M.






If M = 50 and S0 = 15, I0 = 30, (T1)0 = 27 , (T2)0 = 12, N0 = 104 and R0 = 20, then we have the following graphs for the solutions.

Figures 7 and 8 are, respectively, considered for infertile women aged 20-25 and 25-30 with female factor and ovarian reserve of more than 3.5, based on the statistical data of Yazd Infertility enter. In Figure 8, 𝛼 is twice as in Figure 7; therefore, the number of the patients treated by T2 in Figure 8 in its lowest amount (minimum) is half the number of patients treated by T2 in Figure 7. The treatment duration of the patients in Figure 8 is approximately 1.5 times longer than that in patients in Figure 7. Moreover, T1 and T2 in Figure 7, and at the end of the cycle, have an incremental mode. The treatment speed of T2 is more than
that of T1 , since the increase of age, the ovarian reserve decreases, and the medical treatment with clomiphene citrate and gonadotropin will no more be responsive. Therefore, the treatment with IVF T2 is recommended to these patients. The evidences can exactly be observed in Figure 8. Similarly, the recovery speed of R in the Figure 7 reaches its minimum with no longer be responsive, to the extent that the minimum number of the improved patients in Figure 8 is approximately 0.25 of the number of the improved patients in Figure 7. In Figures 7 and 8, when R decreases, the number of patients increases.
Figure 9 has been reviewed for the patients with the ovarian reserve less than 3.5. In Figure 7, u2 is twice more than the u2 in Figure 9. Therefore, the number of the improved patients R in Figure 8 at its lowest value is almost twice as much as than the number of improved individuals in Figure 9. The length of treatment period in Figure 8 is approximately 1.5 times more than the length of treatment in Figure 9, since the speed of referral to the clinic, shown by 𝛾, with respect to the ovarian reserve in Figure 8, is four times more than that in Figure 8.




Figure 7. 𝜆 = 0.47, 𝑏 = 0.34, 𝜂 = 0.47, 𝜇1 = 0.9, 𝑢1 = 0.04, 𝜁 = 0.16, 𝜇2 = 0.92, 𝑢2 = 0.29, 𝛼 = 0.3 and 𝛾 = 0.5.




Figure 8. 𝜆 = 0.47, 𝑏 = 0.34, 𝜂 = 0.47, 𝜇1 = 0.9, 𝑢1 = 0.04, 𝜁 = 0.16, 𝜇2 = 0.92, 𝑢2 = 0.29, 𝛼 = 0.16, and 𝛾 = 0.23.




Figure 9. 𝜆 = 0.47, 𝑏 = 0.34, 𝜂 = 0.47, 𝜇1 = 0.9, 𝑢1 = 0.04, 𝜁 = 0.08, 𝜇2 = 0.96, 𝑢2 = 0.19, 𝛼 = 0.16, and 𝛾 = 0.23.


3. Discussion

The analysis of mathematical models depends on the sustainability of the disease in the community. In this model, the population volume is considered constant at a constant time and provides a good estimate for short periods with similar patients’ conditions. Information obtained from the Yazd Infertility Center confirms that the ovarian reserve decreases as the age increases, thus, the number of patients attending IVF treatment increases, as seen clearly in Figures 7 and 8.

4. Conclusions

In this research, SIT1T2R model is considered to cure infertility in couples. It is proved that the disease free equilibrium point Q0 for SIT1T2R model is locally stable and it is not asymptotically stable, when R0 < 1. Here, R0 is the number of patients who during treatment come across with the illness for the second time, and return to the clinic for secondary treatment. Furthermore, by using the Rung-Kutta method, we solved the model for achieving numerical solution. Considering the numerical results, the patients with ovariarn reserve less than 3.5 with two category of age range 20-25 and 25-30, it is shown that as age
increases the ovarian reserve decreases, hence, Clomiphene Citrate and Gonadotropin treatments are not responsive. Obviously, IVF treatment is recommended in this group of patients. The main achievement of this study was that more ovarian reserve resulted in more cured infertile patients.

Conflict of Interest

The authors declare that they have no competing interests.

Acknowledgments

The authors would like to thank Professor Abbas Aflatoonian for his medical advice.