Mathematical model of sexual orientations in the presence of recovery centers

Introduction

Sexual orientation is the relative attraction one feels toward others (men, women, or both) [1] . Grasping sexual orientation has been a challenging and changing area of research, characterized by notable shifts in societal perceptions and scientific insights over time. Traditionally, sexual orientation was often seen through a binary framework, where heterosexuality was viewed as the standard and other orientations were marginalized or treated as disorders. This view began to shift in the late 1800s and early 1900s, as trailblazers like Magnus Hirschfeld and Alfred Kinsey [2] questioned the dominant beliefs and pushed for a deeper understanding of human sexuality. Since that time, research has increasingly acknowledged sexual orientation as a varied spectrum. Political rights of lesbian, gay, and bisexual (LGB) individuals have significantly improved in numerous Western nations over the past few decades [3] . In the United States, for instance, the Supreme Court determined in June 2015 that states are required to permit marriages between same-sex couples. Nearly 60% of sexual minority groups in the U.S. encountered Sexual Orientation and Gender Identity Change efforts (SOGIECE) from religious environments, counselors, psychologists, and psychiatrists [4] . Mathematical modeling is a method of simulating real-life situations with mathematical equations to predict future behavior by generating a simplified representation of a real system [5, 6] . Such models have been used in evaluating the effects of HIV/AIDS, considering sexual preferences under antiretroviral therapy [7] . Assessing the impacts of bisexual behavior in a global community [8] .

Study the sexual transmission of the Hepatitis B virus (HBV) considering sexual orientation [9] study the spread of HIV/AIDS among LGBT populations [9] . This paper seeks to develop a mathematical model of sexual orientations in the presence of recovery centers.

Model formulation

Classical SIRS (susceptible, infected, recovered and susceptible) is used to model sexual orientations in the presence of recovery centers.

Model assumptions

 It was assumed that homosexuality (same-sex marriage) is not hereditary; rather, people are convinced to practice it.

 The male population can only be influenced by gay and bisexual males to practice homosexuality.  The female population can only be influenced by lesbians and bisexual females to homosexuality.

Model development

The model is formulated into ten compartments; Males (M), Females (F), gays (G), lesbians (L), heterosexual males (Hm), heterosexual females (Hf), bisexual males (Bm), bisexual females (Bf), recovery centers for male (Rm) and recovery centers for females (Rf). Gays, lesbians, bisexual males and bisexual females will be referred to as homosexual populations. Below is the schematic representation of the model.

Invariant regions

The model system 1 will be analyzed in a suitable feasible region as follows.

Boundedness Theorem 1

If the initial conditions of system 1 are in domain Q, then all solutions of system 1 enter and remain in the same domain Q. . Hence the proof.

Positivity

𝑄 = {𝑀(𝑡), 𝐺(𝑡), 𝐿(𝑡), 𝐻 𝑚 (𝑡), 𝐻 𝑓 (𝑡), 𝐵 𝑚 (𝑡), 𝐵 𝑓 (𝑡), 𝑅 𝑚 (𝑡), 𝑅 𝑓 (𝑡)} ∈ 𝑅 + 10 and 𝑀(0) ⩾ 0, 𝐺(0) ⩾ 0, 𝐿(0) ⩾ 0, 𝐻 𝑚 (0) ⩾ 0, 𝐻 𝑓 (0) ⩾ 0, 𝐵 𝑚 (0) ⩾ 0, 𝐵 𝑓 (0) ⩾ 0, 𝑅 𝑚 (0) ⩾ 0, 𝑅 𝑓 (0) ⩾ 0 then, the region will be positively bounded.

Where 𝑪 𝟏 Is a constant of integration. Putting initial conditions

Hence, the region Q is positively invariant, and the model is well-posed mathematically. Since the model monitors human populations, all the variables and parameters are assumed to be non-negative for t ⩾ 0.

Equilibrium states, reproductive number, and stability 4.1 Homosexuality-free equilibrium point (HFE) Theorem 3:

The system 1 has an HFE that is obtained by setting all the homosexual and recovered classes to zero. [11] Proof E 0 = (M 0 ,G 0 ,Hm 0 , Bm 0 , Rm 0 , F 0 , L 0 , Hf 0 , Bf 0 , Rf 0 ) = (ξ1,0, ξ2,0,0, ξ3,0,ξ4,0,0)

The Basic Reproduction Number

The next-generation method [12, [13] was used to obtain reproduction numbers. Using the notation f for a matrix of new infection terms and v for the matrix of the remaining transfer terms in system

It follows that basic reproduction numbers for model 1, denoted by 𝑅 0 𝑚 and 𝑅 0 𝑓 are given by the eigenvalues of F.V -1 as

Local Stability Analysis of Homosexuality Free Equilibrium Point (HFE) Theorem 4:

The HFE of model system 1 is locally asymptotically stable (LAS) if 𝑅 0 𝑚 < 1 and 𝑅 0 𝑓 < 1 and unstable otherwise [12] .

Proof

To establish the local stability, the Jacobian of the model was evaluated at E 0 . The stability of this steady state is then determined based on the signs of the eigenvalues of the corresponding Jacobian, which are functions of the model parameters.

All eigenvalues of (2) were negative under the conditions,

Thus, the homosexuality Free Equilibrium (HFE) of the system is locally asymptotically stable, i.e. 𝑅 0 𝑚 < 1, 𝑅 0 𝑓 < 1 and unstable otherwise.

Global Stability of the homosexuality-free equilibrium point (H.F.E) Theorem 5

If system 1 can be written in the form Where 𝐴 = 𝐷 𝑧 𝐺(𝑋 * , 0) is an M-matrix (the off-diagonal elements of A are non-negative) and Ω is the region where the model makes biological sense [16,17,14] . Then we can conclude that 𝐸 0 is a globally asymptotically stable (GAS) equilibrium of model system 1 provided that 𝑅 0 𝑚 < 1 and 𝑅 0 𝑓 < 1.

Proof

Applying theorem 5 to model system 1 gives

Since 𝐺 ˆ(𝑋, 𝑍) ≥ 0 by theorem 5, 𝐸 0 globally asymptotically stable

Existence of Homosexuality Endemic Equilibrium point (HEEP)

Conditions for the existence of equilibrium are sought for which the model is endemic in the population [12] . For HEEP to exist, the infectious classes of the system 1 must be greater than zero, that is, G > 0, Bm > 0, L > 0, Bf > 0. Setting G, Bm classes to zero and substituting them with each other, and letting them be greater than zero results to The HEEP is globally asymptotically stable whenever 𝑅 0 𝑚 < 1 and 𝑅 0 𝑓 < 1. [15] Proof The following Lyapunov Function was proposed

Where, ∇1 = (κ1 +β1F +π1), ∇2 = (κ1 +ψ1 +µ1F +µ1θ1Bf), ∇3 = (κ1 +γ1G), ∇4 = (κ1 + φ1), ∇5 = (κ1 + η1 + χ1 + Γ1), ∇6 = (κ2 + β2M + π2), ∇7 = (κ2 + ψ2 + µ2M + µ2θ2Bm), ∇8 = (κ2 + γ2), ∇9 = (κ2 + φ2), ∇10 = (κ2 + η2 + χ2 + Γ2).

Letting P * denote all the positive values and Q * denote all the negative values of the solution after substituting evaluated constants, results to Then 𝑑𝑉 𝑑𝑡 = 0 is the singleton 𝐸 * using Lassalle's invariance principle [18, [19] ,

𝑑𝑘(𝑀,𝐺,𝐻 𝑚 ,𝐵 𝑚 ,𝑅 𝑚 ,𝐹,𝐿,𝐻 𝑓 ,𝐵 𝑓 ,𝑅 𝑓 )

Bifurcation Analysis

The Centre Manifold Theorem [16] was used to investigate the nature of the bifurcation of the model 1. Letting 𝑃 𝑘 be the 𝑘 th component of p and

𝑉 𝑘 𝑣 𝑖 𝜕 2 𝑝 𝑘 𝜕𝑦 𝑖 𝛿𝑅 0 (0,0)

Considering a case where 𝑹 𝟎 * 𝒎 = 𝟏 and 𝑹 𝟎 * 𝒇 = 𝟏, 𝒅 𝟏 = 𝒅 𝟏 * and 𝒅 𝟖 = 𝒅 𝟖 * were chosen as bifurcation parameters as follows:

At HFE the Jacobian J * was found to have the right eigenvector given by u =[u1,u2,u3,u4,u5,u6,u7,u8,u9,u10] T, and left eigenvectors as v = [v1, v2, v3, v4, v5, v6, v7, v8, v9, v10]. Solving these vectors resulted to u1 = 0,u2 = 0,u3 = 0,u4 = 0,u5 = 0,u6 = 0,u7 = 0,u8 = 0,u9 =0,u10 = 0 and v1 = 0,v2 = 0,v3 = 0,v4 = 0,v5 = 0, v6 = 0,v7 = 0,v8 = 0,v9 = 0,v10 = 0. Since the values of right and left eigenvectors obtained would give us a = 0 and b = 0, this either means that the system doesn't have a bifurcation, or we can't use the Center Manifold theorem to investigate the bifurcation of this system.

Normalized Sensitivity analysis of the parameters

The technique employed in the sensitivity analysis was based on the reproduction numbers. This helps in determining the level of contribution of each parameter to the reproduction number. This technique uses partial derivatives with respect to the parameters [17] . These indices are given by Λ 𝜌

Partial rank correlation coefficient (PRCC)

To identify parameters that were more influential to the model, PRCC analysis was carried out, where the strength of the outputs and inputs of the model was determined using the LHS method [18] . The closer a value was to +1 or -1, the stronger the linear association between x and y variables, and the farther a value was from +1 or -1, the weaker the linear association between x and y variables.

Numerical Simulations

To observe the dynamics of the sexual orientation model over time, numerical simulations were done using MATLAB software.

The following parameter values were used in the simulation. shows that increasing φ1 (the rate at which the bisexual male population joins recovery centers) by 30% decreases the bisexual male population by 7.9% while increasing φ2 (the rate at which the bisexual female population joins recovery centers) by 30% decreases the bisexual female population by 2.3% at the maximum value. Figure 4 (a) shows that increasing the Bisexual males' recruitment rate (d2) by 69% decreases the heterosexual male population by 17.8%, while increasing the bisexual females' recruitment rate (d5) by 110% decreases the heterosexual female population by 12.3% at maximum value. shows that the population of susceptible (male and female) populations declines with time and that the male population declines faster than the female population.

Results and Discussion

In this study, a deterministic mathematical model of sexual orientation with recovery centers was developed. The qualitative behavior of model 1 was analyzed by obtaining the positivity and boundedness of the solution, equilibrium points, and control reproduction numbers 𝑅 0 𝑚 and 𝑅 0 𝑓 were obtained using the next-generation matrix method. Local stability of HFE was obtained by Routh-Hurwitz criteria, and global stability by the Metzler Matrix method. Lyapunov function was used to determine the global stability of HEEP, and bifurcation analysis was also carried out.

From the analysis carried out, the model was bounded and lay in the positive region. The local stability of HFE was LAS since 𝑅 0 𝑚 < 1 and 𝑅 0 𝑓 < 1, and the global stability of HEEP was asymptotically stable when 𝑅 0 𝑚 > 1 and 𝑅 0 𝑓 > 1. It was also observed that the global stability of HFE was asymptotically stable, given that 𝑅 0 𝑚 < 1 and 𝑅 0 𝑓 < 1. Control reproduction numbers were obtained as 𝑅 0 𝑚 = 3.9 and 𝑅 0 𝑓 = 3.2. Moreover, a sensitivity analysis of the system was also performed based on the control reproduction number, and it was observed that d1, d8 were the most influential parameters hence the conclusion that varying d1,d8 would significantly influence the spread of homosexuality within the community.

To identify the highly sensitive parameters in the model, the sensitivity analysis was performed based on the PRCC method and LHS, and showed (e1,e3, Γ2,χ1,χ2) as the most influential parameters. Finally, numerical simulations showed that increasing recruitment rates (d1,d2,d5,d8) increased homosexuality in the population and increasing 𝜑 1 and 𝜑 2 decreased bisexual male and bisexual female population, respectively. Therefore, this paper shows that the best way of reducing the spread is by reducing the contact rate and increasing transfer rates from homosexual classes to recovery centers.

Acknowledgement

AcknowledgmentThe authors appreciate the support accorded by the Kirinyaga University staff fraternity.

References

1. Bailey JM, Vasey PL, Diamond LM, Breedlove SM, Vilain E, Epprecht ME. Sexual orientation, controversy, and science. Psychological Science in the Public Interest. 2016;17(2):45-101.
2. Pillard RC, Bailey JM. Human sexual orientation has a heritable component. Hum Biol. 1998;70(3):347-365.
3. Carroll A, Itaborahy L. State sponsored homophobia 2015: a world survey of laws: criminalisation, protection and recognition of same-sex love. Geneva: International Lesbian, Gay, Bisexual, Trans and Intersex Association; 2015.
4. Turban JL, King D, Reisner SL, Keuroghlian AS. Psychological attempts to change a person's gender identity from transgender to cisgender: estimated prevalence across US States, 2015. Am J Public Health. 2019;109(10):1452-1454.
5. Abuasbeh K, Qureshi S, Soomro A, Awadalla M. An optimal family of block techniques to solve models of infectious diseases: Fixed and adaptive stepsize strategies. Mathematics. 2023;11(5):1135.
6. Ramos H, Qureshi S, Soomro A. Adaptive step-size approach for Simpson's-type block methods with time efficiency and order stars. Computational and Applied Mathematics. 2021;40(6):219.
7. Sriningsih R, Putra AA, Yuniarti E, Soleh M, et al. Construction of mathematical model between HIV-AIDS and lesbian, gay, bisexual, and transgender (LGBT) transmission in a population. In: Journal of Physics: Conference Series, 2020, p. 12055.
8. Espitia CC, Botina MA, Solarte MA, Hernandez I, Riascos RA, Meyer JF. Mathematical model of HIV/AIDS considering sexual preferences under antiretroviral therapy, a case study in San Juan de Pasto, Colombia. J Comput Biol. 2022;29(5):483-493.
9. Lu M, Shu Y, Huang J, Ruan S, Zhang X, Zou L. Modelling homosexual and heterosexual transmissions of hepatitis B virus in China. J Biol Dyn. 2021;15(1):177194.
10. Birkoff G, Rota GC. Ordinary differential equations. Boston: Ginn, 1982.
11. Ngari C, Yrusgitonga S. Modelling the impact of screening, treatment and underlying health conditions on dynamics of COVID-19. 2023 [cited 2024].
12. Yang HM. The basic reproduction number obtained from Jacobian and next generation matrices-A case study of dengue transmission modelling. Biosystems. 2014;126:52-75.
13. Castillo-Chavez C, Feng Z, Huang W. On the computation of R0 and its role on global stability, 2001.
14. Khalil KH. Lyapunov stability. In: Control systems, robotics and automation. 2009; p. 115.
15. Hassard B, Wan YH. Bifurcation formulae derived from center manifold theory. J Math Anal Appl. 1978;63(1):297-312.
16. Petzold L, Li S, Cao Y, Serban R. Sensitivity analysis of differential-algebraic equations and partial differential equations. Comput Chem Eng. 2006;30(10-12):15531559.
17. Gasior KI. Examining the influence of nondimensionalization on partial rank correlation coefficient results when modeling the epithelial mesenchymal transition. Bull Math Biol. 2025;87(1):153. DOI: 10.1007/S11538-024-01393-Y/FIGURES/37
18. Oluwafemi T, Bako D. Numerical solution of a mathematical model of the effect of same-sex marriage. 2022;1:5-7.