Qualitative analysis of a mathematical model of divorce epidemic with anti-divorce therapy

1. Introduction

Marriage is a socially accepted union or legal contract between two persons that established rights and obligations between them and their in-laws as well as the society. It can sometimes be referred to as the union of matrimony or wedlock, where interpersonal sexual relationship is biologically acknowledged [1]. The institution of marriage originally was meant to be a union of no bitterness. However, ugly experiences have emerged in which divorce and separation remain the most re-occurring decimals. The later phenomena are likely to happen when the ratio of negative to positive behaviour is equal or greater than unity. Divorce is the dissolution of marriage between two persons [2] which has a deformable mark on the family structure in our present era. This bitter trend continues to persist across Africa and the developed nations with no stigmatization. The immediate causes of divorce may not be limited to early marriage, less education, low income, premarital cohabitation and pregnancy, infidelity and lack of religion affiliation [3].

Divorce is an endemic issue that seriously affects the social and economic structure of contemporary society as much as any disease. Approximately one – half of all first marriages in the world end up in divorce or separation [4], with even higher rates of divorce for second marriages [5]. Separation leads to divorce as 75% of separation eventually result in divorce [6]. Divorce terms to increase childhood poverty and illiteracy rates [7]. It becomes important to eliminate this deadly social virus from our society.

Mathematical modelling has been used as a veritable tool for the control of epidemics in the past decades, and we can adopt it for the containment of the spread of divorce in marital institutions. In the past, concern has been focused on the influence of economy [8] and social contagion [9] on divorce models. Recently, [2] formulated a model for the spread of divorce in Ghana with three major compartments; Married, Divorced and Separated cases. Meanwhile, [10] extended their work to include the population of singles without any divorce prevention measure. Other mathematical studies on this subject matter can be found in [11,12,13]. Like other authors [2,9,10], we constructed a mathematical model for the spread of divorce epidemic. Apart from that, we incorporate a class of restored marital cases with anti-divorce therapy and reconciliation to the model of [10] using standard incidence function as the force of marital disunity, where the social family disorder is spread by divorced and separated persons over the married individuals. This is a major feature missing in the past works.

The organization of this paper is as follows: In Section 2, we formulate the model and establish the existence of equilibra and their local stability. In Section 3, we prove the global dynamics of each of the feasible steady states of the system (2) by constructing Lyapunov functions. In Section 4, we present the results and discussion of the numerical simulation. Finally, a brief conclusion is given in Section 5 to end this work.

2. Model construction

The present model derived its motivation from the work of [10], in which for the purpose of clarity, the basic model can be stated as; A more realistic model of divorce epidemic (1) with anti-divorce therapy is formulated based on the assumptions below and the schematic diagram in Figure 1.

1. People who are single or individuals who are ready for marriage but not yet married.
2. Some Separated cases may resort to divorced but not vice-versa.
3. Divorce or Separated cases restored can remain unbroken.
4. Only Separated cannot be remarried.
5. People who divorce can remarry or remain single.
6. The anti-divorce parameter changes from 0 to 1 \((0\leq\phi\leq1)\).

Table 1. Variables and parameters of the model (2).

<b>Variables</b>	<b>Description</b>
\(S_{1}(t)\)	Number of singles who are due for marriage at time \(t\)
\(M(t)\)	Number of Married cases at time \(t\)
\(D(t)\)	Number of Divorced cases at time \(t\)
\(S_{2}(t)\)	Number of Separated cases  at time \(t\)
\(R(t)\)	Number of Restored Marital cases at time \(t\)
Parameters	Description
\(Q\)	Recruitment rate for the Single individuals
\(\beta_{1}\)	Divorced rate of the married
\(\beta_{2}\)	Separated rate of the married
\(\alpha_{1}\)	Rate of getting married by the singles
\(\alpha_{2}\)	Rate of re-marriage after divorced
\(d_{3}\)	Separated rate of the married
\(d_{2}\)	Proportion of married individuals that divorced
\(d_{1}\)	Proportion of married cases that results into Separation
\(r_{1}\)	Rate of restoring divorced cases in marriage by reconciliation
\(r_{2}\)	Rate of restoring Separated cases in marriage through reconciliation
\(\gamma\)	Rate of remaining single after divorce
\(\phi\)	Anti-divorce parameter \((0\leq\phi\leq1)\)
\(\mu\)	Natural death rate of individuals
\(\lambda_{m}\)	Force of marital disunity

Note that, the first four Variables and the first seven Parameters in Table 1 have the same meaning as in model (1) [10]. However, \(\beta\) which is the rate of Separated getting married again in model (1) is disputed and so ignored in our study.

From the model flow-diagram of divorce epidemic and assumptions, the following differential equations are derived;

where Equation (3) is the force of marital disunity. The sum of the entire system (2) yield with \(N=S_{1} +M+D+S_{2}+R.\) The system (2) can be studied within the feasible domain \begin{equation*} D_{m}=\left\lbrace S_{1}>0, M\geqslant0, D\geqslant0, S_{2}\geqslant0, R\geqslant0\mid N=\frac{Q}{\mu}\right\rbrace. \end{equation*}

2.1. Model steady states

At the steady state, the system (2) takes the form where \(\omega_{1}=\alpha_{1} +\mu, \omega_{2}=\mu+\gamma+r_{1}+\alpha_{2}\) and \(\omega_{3}=\mu+r_{2}+d_{3}.\) The solution of the system (5) gives the endemic equilibra in terms of \( \lambda_{m}^{**}\). with Fundamentally, the divorce-free equilibrium (DFE) can be obtained from the relation (7) above. Thus, introducing the expressions (6) into (7) and simplifying accordingly we arrive at such that Clearly from (8) \(\lambda_{m}^{**}=0 \), which is significant for evaluating DFE. Consequently substituting \(\lambda_{m}^{**}=0 \) into the expressions (6) give \begin{equation*} E^{0}=(S_{1}^{0}, M^{0}, D^{0}, S_{2}^{0}, R^{0})=\left(\frac{Q}{\alpha_{1}+\mu},\frac{Q \alpha_{1}}{\alpha_{1}+\mu}, 0, 0, 0\right). \end{equation*}

2.2. Divorce control reproduction number

The computation of the divorce controlled reproduction number denoted here by \(R_{ed} \) is done following the standard next generation matrix explained in [14] with the usual notation being given as; \[R_{ed}=\rho \left(F V^{-1}\right)=(1-\phi)\frac{\alpha_{1}}{\omega_{1}} \left(\frac{\beta_{1} (d_{1} \omega_{3}+ d_{2} d_{3})+\beta_{2} d_{2} \omega_{2}}{\omega_{2} \omega_{3}}\right),\] with \[F=(1-\phi)\frac{\alpha_{1}}{\omega_{1}}\begin{pmatrix} d_{1}\beta_{1} & d_{1}\beta_{2} \\ d_{2}\beta_{1} & d_{2}\beta_{2}\\ \end{pmatrix}\quad \text{and}\quad V^{-1}=\begin{pmatrix} \frac{1}{\omega_{2}} & \frac{d_{3}}{\omega_{2} \omega_{3}}\\ 0 & \frac{1}{\omega_{3}}\\ \end{pmatrix}.\] In the absence of anti-divorce therapy and reconciliation, \( \phi=0, r_{1}=r_{2}=0 \), we derive divorce reproduction number as; \[R_{d}=\frac{\alpha_{1}}{\omega_{1}} \left[\frac{\beta_{1} [d_{1} (\mu+d_{3})+d_{2} d_{3}]+\beta_{2} d_{2} (\mu+\gamma+\alpha_{2})}{(\mu+d_{3})(\mu+\gamma+\alpha_{2})}\right].\] However, for single intervention in terms of anti-divorce therapy and reconciliation only we have the respective reproduction numbers as; \begin{align*} R_{ad}&=(1-\phi)\frac{\alpha_{1}}{\omega_{1}} \left[\frac{\beta_{1}[d_{1} (\mu+d_{3})+d_{2} d_{3}]+\beta_{2} d_{2} (\mu+\gamma+\alpha_{2})}{(\mu+d_{3})(\mu+\gamma+\alpha_{2})}\right],\\ R_{ec}&=\frac{\alpha_{1}}{\omega_{1}} \left(\frac{\beta_{1} (d_{1} \omega_{3}+ d_{2} d_{3})+\beta_{2} d_{2} \omega_{2}}{\omega_{2} \omega_{3}}\right).\end{align*} Interestingly, it can be deduced from Theorem 2 in [14] that:

Claim 1. The DFE is locally asymptotically stable in the case \( R_{ed}< 1\) and unstable when \( R_{ed}>1\).

2.3. Existence and local stability of divorce-endemic equilibrium

The existence of divorce endemic equilibrium follows immediately from (8) where \(\lambda_{m}^{**}\neq0\). Therefore, solving the remaining part of the equation leads to Thus, the divorce-endemic equilibrium denoted by DE after substituting (10) into (6) is given by \(DE=(S_{1}^{**},M^{**},D^{**},S_{2}^{**},R^{**} )\).

Recall from (10) that \(\lambda_{m}^{**}>0\) iff \(R_{ed}>1\) and \(\alpha_{1}>\alpha_{2}\), which are the necessary conditions for divorce to persist in marriage institution.

The method of linearization process of the system (2) at DE gives the following Jacobian Matrix

From the Jacobian matrix in (11), \( -\mu \) is an eigenvalue and the remaining ones are gotten from the matrix below with the corresponding characteristics equation defined as follows where \begin{align*} a_{4} =&1,\\ a_{3} =& \mu+\left(\sum_{i=1}^{3}\omega_{i}\right)+\lambda_{m}^{**}-(1-\phi) (d_{1} \beta_{1}+d_{2} \beta_{2})\left(\frac{M^{**}}{N^{**}}\right),\\ a_{2}=&\lambda_{m}^{**}\left[\left(\sum_{i=1}^{3}\omega_{i}\right)-d_{1} \alpha_{2}\right]+\mu\left(\sum_{i=1}^{3}\omega_{i}\right)+\omega_{1}\left(\sum_{i=2}^{3}\omega_{i}\right)-d_{2} (1-\phi)^{2} \left(\frac{M^{**}}{N^{**}}\right)^{2} \beta_{1}\beta_{2}(1-d_{1})\\ &+\omega_{2} \omega_{3}-(1-\phi)\left(\frac{M^{**}}{N^{**}}\right)\left[\beta_{2} d_{2} \omega_{2} +\beta_{!} (d_{1}\omega_{3}+d_{2} d_{3})+(\mu+\omega_{1})(d_{1}\beta_{1}+d_{2} \beta_{2})\right],\\ a_{1} =&\omega_{1}\omega_{2}\omega_{3} +\mu\left[\omega_{1}\left(\sum_{i=2}^{3}\omega_{i}\right)+\omega_{2}\omega_{3}\right]+\lambda_{m}^{**}\left[\omega_{1}\left(\sum_{i=2}^{3}\omega_{i}\right)-d_{2}\alpha_{2}(1-\phi) \left(\frac{M^{**}}{N^{**}}\right)\right.\beta_{2}(1-d_{1})\\ &+\omega_{2}\omega_{3}-d_{1}(\gamma\alpha_{1}+\alpha_{2}\omega_{1})-\alpha_{2}(d_{1}\omega_{3}+d_{2} d_{3})\Bigg]+(1-\phi)^{2} \left(\frac{M^{**}}{N^{**}}\right)^{2} \beta_{1}\beta_{2}\left[2\mu+\alpha_{2}(1-d_{1})\right]\\ &-(1-\phi)\left(\frac{M^{**}}{N^{**}}\right)\left[\beta_{1}(\alpha_{1}d_{2}d_{3}+d_{1}(\omega_{1}\omega_{3}+\mu(\omega_{1}+\omega_{3})))+\beta_{2}d_{2}(\mu\left(\sum_{i=1}^{2}\omega_{i}\right)+\omega_{2}\omega_{3})\right],\\ a_{0} =&\lambda_{m}^{**}\left[\omega_{1}\omega_{2}\omega_{3}-(\gamma\alpha_{1}+\alpha_{2}\omega_{1})d_{2}(1-\phi)\left(\frac{M^{**}}{N^{**}}\right)\beta_{2}(1+d_{1})+(d_{1}\omega_{3}+d_{2}d_{3})\right]\\ &+\mu\omega_{1}\left[\omega_{2}\omega_{3}-(1-\phi)\left(\frac{M^{**}}{N^{**}}\right)(\beta_{2}d_{2}\omega_{2}+\beta_{1}(d_{1}\omega_{3}+d_{2}d_{3}))-d_{2}(1-\phi)^{2}\left(\frac{M^{**}}{N^{**}}\right)\beta_{1}\beta_{2}(1-d_{1})\right]. \end{align*} To guarantee that all roots of (13) are real and negative, the Routh-Hurwitz stability criterion [15] requires; Thus, it is evident to see that all the coefficient of (13) are strictly positive if \[\lambda_{m}^{**}>0\;\; (\Leftrightarrow R_{ed}>1)\ \ \text{and}\ \ \frac{M^{**}}{N^{**}} < 1.\] Therefore, the inequalities (14) hold, and we cam claim that:

Claim 2. The divorce endemic equilibrium DE is locally asymptotically if \( R_{ed}>1\) and \(\frac{M^{**}}{N^{**}}< 1\).

3. Global dynamics

3.1. Global stability of DFE

This section discusses the global behavior of the model (2) at the divorce free equilibrium state following the well-known stability method of Lyapunov functions [16].

Claim 3. The system (2) admits a globally asymptotically DFE iff \( R_{ed}< 1\).

Proof. By matrix theoretic approach as done in [16] we derive the following Lyapunov function \[L(t)=\frac{\beta_{1}}{\omega_{2}} D +\left(\frac{\beta_{i} d_{3}}{\omega_{2} \omega_{3}}+\frac{\beta_{2}}{\omega_{3}}\right) S_{2},\] whose time derivative is given by \begin{align*} L^{'}(t) =&\frac{\beta_{1}}{\omega_{2}}\left[d_{1} \left((1-\phi) \frac{\beta_{1} D + \beta_{2} S_{2}}{N}\right)M + d_{3} S_{2}-\omega_{2} D\right]\\ &+\left(\frac{\beta_{i} d_{3}}{\omega_{2} \omega_{3}}+\frac{\beta_{2}}{\omega_{3}}\right)\left[d_{2}\left((1-\phi) \frac{\beta_{1} D + \beta_{2} S_{2}}{N}\right)M -\omega_{3} S_{2}\right]. \end{align*} Collecting like terms and simplifying leads to \begin{align*} L^{'}(t) & =(\beta_{1} D + \beta_{2} S_{2})\left[(1-\phi) \left(\frac{\beta_{1} (d_{1} \omega_{3}+ d_{2} d_{3})+\beta_{2} d_{2} \omega_{2}}{\omega_{2} \omega_{3}}\right)\frac{M}{N}-1\right]. \end{align*} Recall that at divorce-free equilibrium we have \(\frac{M}{N}=\frac{M^{0}}{N^{0}}=\frac{\alpha_{1}}{\omega_{1}}\) and so \begin{align*} L^{'}(t) & =(\beta_{1} D + \beta_{2} S_{2})\left[(1-\phi)\frac{\alpha_{1}}{\omega_{1}} \left(\frac{\beta_{1} (d_{1} \omega_{3}+ d_{2} d_{3})+\beta_{2} d_{2} \omega_{2}}{\omega_{2} \omega_{3}}\right)-1\right]. \end{align*} Thus, \( L^{'}(t) =(\beta_{1} D + \beta_{2} S_{2})(R_{ed}-1), \), \(L^{'}(t)\leqslant0\) if \(R_{ed}\leqslant1\) and \( D=S_{2}=0\). Therefore, by Lyapunov's stability theory [15], the equilibrium point \(E^{0} \) is globally asymptotically stable.

3.2. Global stability of DE

Claim 4. The divorce-endemic equilibrium of the system (2) at \(\alpha_{2}=0\) is globally stable if \(R_{ed}>1\).

Proof. Taking the Lyapunov function \( W_{1}(t)\) as in [17]

with \( m_{1}, m_{2}, . . . , m_{5} \) are positive constants to be properly chosen.

The differential coefficient of \(W_{1} (t)\) in (15) with the substitution of model (2) at \( \alpha_{2}=0\) gives

At divorce endemic equilibrium, we have \begin{align*} &\left(\alpha_{1} +\mu\right)S_{1}^{**} =Q+\gamma D^{**},\\ &\alpha_{1}S_{1}^{**}= \mu M^{**}+\left(\frac{\beta_{1} D^{**} + \beta_{2} S_{2}^{**}}{N}\right)M^{**},\\ &\left(\mu+\gamma+r_{1}\right)D^{**}=d_{3} S_{2}^{**}+d_{1} \left(1-\phi\right) \left(\frac{\beta_{1} D^{**} + \beta_{2} S_{2}^{**}}{N}\right)M^{**},\\ &\left(\mu+r_{2}+d_{3}\right)S_{2}^{**}=d_{2}\left(1-\phi\right) \left(\frac{\beta_{1} D^{**} + \beta_{2} S_{2}^{**}}{N}\right)M^{**},\\ &\mu R^{**} =r_{1} D^{**}+ r_{2} S_{2}^{**}. \end{align*} Therefore, where \begin{align*} G \left(S_{1}, M, D, S_{2}, R\right) =&-m_{1}\left(1-\frac{S_{1}^{**}}{S_{1}}\right)\left(1-\frac{D}{D^{**}}\right)\gamma D^{**}+m_{2}\left(1-\frac{M^{**}}{M}\right)\left(\frac{S_{1}}{S_{1}^{**}}-\frac{M}{M^{**}}\right)\mu M^{**}\\ & +m_{2}\left(1-\phi\right)\beta_{1}\left(1-\frac{M^{**}}{M}\right) \left(\frac{S_{1}}{S_{1}^{**}}-\frac{D M}{D^{**} M^{**}}\right)\frac{D^{**} M^{**}}{N}\\ & +m_{2}\left(1-\phi\right)\beta_{2}\left(1-\frac{M^{**}}{M}\right) \left(\frac{S_{1}}{S_{1}^{**}}-\frac{S_{2} M}{S_{1}^{**} M^{**}}\right)\frac{S_{2}^{**} M^{**}}{N}\\ & +m_{3} d_{1}\left(1-\phi\right)\beta_{1}\left(1-\frac{D^{**}}{D}\right) \left(\frac{D M}{D^{**} M^{**}}-\frac{D}{D^{**}}\right)\frac{D^{**} M^{**}}{N}\\ & +m_{3} d_{1}\left(1-\phi\right)\beta_{2}\left(1-\frac{D^{**}}{D}\right) \left(\frac{S_{2} M}{S_{2}^{**} M^{**}}-\frac{D}{D^{**}}\right)\frac{S_{2}^{**} M^{**}}{N}\\ & +m_{3} d_{1}\left(1-\phi\right)\beta_{2}\left(1-\frac{D^{**}}{D}\right) \left(\frac{S_{2}}{S_{2}^{**}}-\frac{D}{D^{**}}\right)d_{3} S_{2}^{**}\\ & +m_{4} d_{2}\left(1-\phi\right)\beta_{1}\left(1-\frac{S_{2}^{**}}{S_{2}}\right) \left(\frac{D M}{D^{**} M^{**}}-\frac{S_{2}}{S_{2}^{**}}\right)\frac{D^{**} M^{**}}{N}\\ & +m_{4} d_{2}\left(1-\phi\right)\beta_{2}\left(1-\frac{S_{2}^{**}}{S_{2}}\right) \left(\frac{S_{2} M}{S_{2}^{**} M^{**}}-\frac{S_{2}}{S_{2}^{**}}\right)\frac{S_{2}^{**} M^{**}}{N}\\ & +m_{5}\left(1-\frac{R^{**}}{R}\right) \left(\frac{D}{D^{**}}-\frac{R}{R^{**}}\right) r_{1} D^{**} +m_{5}\left(1-\frac{R^{**}}{R}\right) \left(\frac{S_{2}}{S_{2}^{**}}-\frac{R}{R^{**}}\right) r_{2} S_{2}^{**}, \end{align*} and \(G\) is a non-positive expression as supported by the Lemmas 2.2, 2.3 [18] and [17].

Since \((S_{1},M, D, S_{2},R)\neq(S_{1}^{**},M^{**},D^{**},S_{2}^{**},R^{**})\), \(W_{1}^{'} (t)< 0\), \(W_{1}^{'} (t)=0\) if \(2-\frac{S_{2}}{S_{2}^{**}}=\frac{S_{2}^{**}}{S_{2}}\leqslant0\) and \(S_{1}=S_{1}^{**}, M=M^{**}, D=D^{**}, S_{2}=S_{2}^{**}, R=R^{**} \). Hence, DE is globally stable by Lyapunov's stability theory [15].

4. Results and discussions

To aid the understanding of the analyzed results, we performed the numerical simulation of the proposed model using MATLAB ode45 solver and the parameter values given in Table 2.

Table 2. Values of variables and parameters of the model (2) used for numerical simulations.

<b>Initial Variables</b>	<b>Value</b>	<b>Source</b>
\(S_{1}(0)\)	50	Assumed
\(M(0)\)	40	Assumed
\(D(0)\)	15	Assumed
\(S_{2}(0)\)	15	Assumed
\(R(0)\)	10	Assumed
<b>Parameter</b>	<b>Value</b>	<b>Source</b>
\(Q\)	2	Assumed
\(\beta_{1}\)	0.022	\cite{j}
\(\beta_{2}\)	0.031	\cite{j}
\(\alpha_{1}\)	0.101	\cite{j}
\(\alpha_{2}\)	0.061	\cite{j}
\(d_{3}\)	0.021	\cite{j}
\(d_{2}\)	0.7	Assumed
\(d_{1}\)	0.3	Assumed
\(r_{1}\)	0.1	Assumed
\(r_{2}\)	0.3	Assumed
\(\gamma\)	0.01	Assumed
\(\phi\)	0.2	Assumed
\(\mu\)	0.02	Estimated

In Figure 2(a), we observe that married couples develop the potentials to stay or remain in marriage longer when they attained marriage seminars/courses and exhibit the spirit of reconciliation in handling their differences. On the other hand, the union can easily break apart when the above mentioned controls are ignored.

The impact of anti-divorce control is depicted on the divorced cases in Figure 2(b). Here, the cases of marriage breaking apart (divorce) persist uniformly in the family circle where reconciliation and anti-divorce therapy are absent. Furthermore, the results supports that with the controls, divorced cases can be eliminated in marriage institutions and marriage couples who stay in marriage up to 30 years may have a stable family structure and no longer divorced their partners. This agrees with the result of high commission for planning of Morocco on marriage and divorce as contained in [11] which says that, couples with 20 years of experiences has a decreased divorce rate at 3%.

In the case of separated in Figure 2(c), the number rises higher for the scenario where the couples have not had proper premarital education and do not reconcile their misunderstanding on continuous bases. More so, the study reveals that marriage couples that live for over 15 years without separation are less venerable to divorce tendencies in marriage.

Figure 3(a) is a simple demonstration that increasing the rate of reconciliation impacts positively on the broken families by reuniting and hence increasing the number of restored cases. By implication, separated or divorced families will remain broken without adequate reconciliation process. Meanwhile, the effect of \( \alpha_{2}\) on the divorced is illustrated in Figure 3(b). The positive variation in \( \alpha_{2}\) reduces the cases of divorce in marriage. Therefore, it is recommendable to say that the re-marriage among the divorced helps to re-unite broken families.

The separated/divorced cases grow exponentially with respect to the married cases in Figure 4(a). This means that as long as married couples lived together, marital disorder in terms of separation or divorce may always be experienced among them. In other words, divorce cannot occur without marriage and the more marriage cases, the large scale of divorced population expected.

In Figure 4(b) we observed the behavioral trends of the model reproduction numbers that follows the inequality \(R_{ed}< R_{ec}< R_{ad}< R_{d}\). This symbolizes that the combination of reconciliation and anti-divorced protocols is beneficial in stabilizing marriages than when a single intervention is used. Furthermore, it is important to emphasize that, reconciliation in family settlement is better than marriage seminar/courses as proposed in [10]. This is also in line with the result of [19] that says; therapy is only a support but not a major component in solving marriage problems since \(R_{ec}< R_{ad} \). However, \(R_{d}\) standing above the rest of the reproduction ratios means that marriages without controls are fragile to crack or break apart in future.

5. Conclusion

A non-linear ordinary differential system of equations for examining the spread of divorce epidemic with anti-divorce therapy has been proposed and analyzed. Important results from the qualitative analysis of the modified model reveal that, the model has globally stable equilibrium states, namely the divorce-free and endemic equilibrium. Numerical results of this study suggest that the presence of correctional practices such as marriage seminar/courses and reconciliation efforts can ensure prolonged stable marriages and prevents cracks and repairs the breaking points (separation or divorce) in family structure. In addition, the marriage couples that escape separations in early 15 years of marriage and stay together in the next 15 years cannot divorce any longer.

Author Contribution: All authors contributed equally to the writing of this paper. All authors read and approved the final manuscript.

Conflicts of Interest: The authors declare no conflict of interest.