In this paper, we derive and analyze a mathematical model of a sexual response. As a starting point, we discuss two studies that proposed a connection between a sexual response cycle and a cusp catastrophe and explain why that connection is incorrect but suggests an analogy with excitable systems. This then serves as a basis for derivation of a phenomenological mathematical model of a sexual response, in which the variables represent levels of physiological and psychological arousal. Bifurcation analysis is performed to identify stability properties of the model’s steady state, and numerical simulations are performed to illustrate different types of behavior that can be observed in the model. Solutions corresponding to the dynamics associated with the Masters–Johnson sexual response cycle are represented by “canard”-like trajectories that follow an unstable slow manifold before making a large excursion in the phase space. We also consider a stochastic version of the model, for which spectrum, variance, and coherence of stochastic oscillations around a deterministically stable steady state are found analytically, and confidence regions are computed. Large deviation theory is used to explore the possibility of stochastic escape from the neighborhood of the deterministically stable steady state, and the methods of an action plot and quasi-potential are employed to compute most probable escape paths. We discuss implications of the results for facilitating better quantitative understanding of the dynamics of a human sexual response and for improving clinical practice.

## I. INTRODUCTION

Over the last several decades, mathematical models have successfully been used to study a variety of problems in human physiology, from respiration to heart and blood circulation, endocrine, renal and gastrointestinal systems, hearing and inner ear, retina and vision, and muscles.^{1–3} At the same time, with sex being an essential part of human existence since the dawn of civilization, mathematical modeling of a sexual response has so far been almost non-existent. This can largely be attributed to sex being a taboo subject, which for centuries hindered collection of comprehensive and accurate data on sexual practices and the dynamics of a sexual response.

Significant progress in this direction was made in the twentieth century, starting with the works of Freud on origins of sexuality^{4} and continuing with an extensive work by Kinsey and his colleagues on collecting and producing comprehensive surveys of diversity and frequency of various human sexual practices.^{5,6} Despite subsequent criticisms of some of their methodology, in many respects, Kinsey reports have stood the test of time and are still regularly used as sources of data on various aspects of human sexuality. The next major breakthrough was made by Masters and Johnson,^{7} who introduced the notion of a *human sexual response cycle*, a sequence of physiological changes taking place in the organism during sex. They proposed a four-stage model of a sexual response cycle: excitation $\u2192$ plateau $\u2192$ orgasm $\u2192$ resolution (EPOR model) and described detailed physiological changes associated with each stage of the cycle in males and females.

As already noted by Masters and Johnson, there are significant differences between male and female sexual response cycles, including, among others, the presence of a refractory stage for males, during which additional stimulation cannot produce new arousal, as well as the possibility of multiple orgasms in females. To account for these differences, a number of alternative models of a sexual response cycle were subsequently proposed, including a desire–excitation–orgasm–resolution (DEOR) model,^{8–11} a “sexual man model,”^{12} and the dual control model.^{13} Some of the criticisms of the original Masters–Johnson raised in these later works included lack of contribution from the psychological component of the arousal^{14–16} and an observation that arousal was continuing to increase during the plateau phase.^{17} Levin^{18,19} provides a nice overview and comparison of these and other models of a human sexual response.

Substantial progress has been made in the last 20 years in terms of measuring various physiological reactions during a sexual response cycle using fMRI.^{20–22} Besides providing valuable data on temporal physiological changes, including activation of different parts of the brain during sexual stimulation and orgasms, these studies have further highlighted differences in the sexual response of males and females.

In this paper, we propose and study a phenomenological mathematical model of a sexual response cycle aimed at reproducing and explaining salient features of the Masters–Johnson EPOR model in human males. To overcome some of the above-mentioned deficiencies of that model, we will represent sexual arousal as a system of two coupled ordinary differential equations, with the two variables representing current levels of physiological and psychological arousal during sex. To provide some further background to the model, in Sec. II, we review two series of studies that pointed to a tenuous connection between a sexual response cycle and catastrophe theory, and we explain why that characterization is not correct. Then, in Sec. III, we present a derivation of the mathematical model of a sexual response, for which stability analysis and numerical simulations are performed in Sec. IV. Section V considers the role of stochastic effects, including analysis of stochastic fluctuations around a deterministically stable steady state, and the calculation of most probable stochastic escape paths. This paper concludes in Sec. VI with a discussion of results.

## II. AN ORGASM IS *NOT* A CATASTROPHE

^{23,24}to evolutionary theory

^{25,26}and sports performance.

^{27,28}Two series of papers, one by Hubey

^{29,30}and one by Levin,

^{18,31}have suggested that the human sexual response curve and, in particular, orgasms are associated with the so-called

*cusp catastrophe*. To explain why in both cases this characterization is incorrect, first, we briefly review cusp catastrophe as proposed by Zeeman.

^{32,33}Consider a one-dimensional gradient dynamical system

*cusp catastrophe*: as the values of $b$ increases and crosses zero, there is a cusp in the $a$– $b$ parameter plane, with the system having one stable steady state before the cusp and two stable equilibria after the cusp.

^{29,30}suggested that a sexual response can be mathematically modeled as a Duffing oscillator,

^{34}

Hubey^{29,30} plotted an amplitude response curve for the model (2), i.e., amplitude of oscillations $ |\theta |$ around zero steady state of the unforced model, as a function of forcing frequency $\Omega $, and then drew two conclusions: a sexual response cycle is a periodic trajectory, and there is a cusp catastrophe in the model due to the presence of sudden jumps between upper and lower branches on the amplitude response curve. The first of these conclusions is incorrectly associating motion along the amplitude response curve with periodic solutions represented as closed curves in the phase plane for two-dimensional systems. In the former case, once the forcing frequency is fixed, this will result in some fixed amplitude that will not be changing with time. To move along the amplitude response curve, one would actually have to change a forcing frequency, and then for each value of the frequency, there would be an associated value of the amplitude of oscillations. In contrast, closed phase trajectories in the phase plane do represent periodic solutions once system parameters permit existence of such solutions for certain parameter values. There is no cusp catastrophe associated with the tip of the amplitude response curve, because even though at that point, a transition between branches occur, in each case, the system would exhibit oscillations with only one value of the amplitude, and there is no change in the number of steady states of the model. Moreover, whereas the Duffing oscillator can exhibit periodic, quasiperiodic, and chaotic behavior, as shown in Fig. 2, *none of its solutions* has the form that would qualitatively resemble temporal dynamics of the Masters–Johnson sexual response cycle. Furthermore, the Masters–Johnson sexual response cycle is concerned with a single trajectory showing growth of arousal culminating in an orgasm, while any point on the amplitude response curve rather corresponds to sustained periodic oscillations of a certain fixed amplitude. All this suggests that while there may be some scope in modeling the dynamics of sexual arousal as a forced two-dimensional system, associating hysteresis on the amplitude response curve with a cusp catastrophe is not justified.

Levin^{18,31} has suggested a modification of a Masters–Johnson graphical representation of a sexual response cycle, which addresses the deficiency of the original figure in terms of correctly identifying the refractory period as the period of time after an orgasm, during which subsequent stimulation cannot lead to a new arousal. Additionally, it also introduced a new feature; namely, at the point of an orgasm, the arousal is assumed to instantaneously change its direction from growing to decreasing, and this was the argument used to suggest the existence of a cusp catastrophe. There are two reasons why this interpretation is incorrect. The first one concerns the fact that the cusp catastrophe is associated with a change in the number of steady states of the model, not with changes in its dynamical behavior, however rapid they may be. The second issue is that looking carefully at this figure, we note that the derivative of central arousal is discontinuous at the time point of an orgasm, whereas dynamical systems, for which the cusp catastrophe is defined, are described by differentiable vector fields, for which such discontinuity would be impossible. Our interpretation of this figure is that it is rather a conceptual representation of separation of time scales, with a slow time scale during the plateau phase, followed by a fast time scale of resolution after an orgasm, which resembles the dynamics exhibited by excitable systems, and this is exactly the approach we will use in our model presented in Sec. III.

## III. MODEL DERIVATION

^{35}that a sexual response is a complex process with multiple contributions, we model sexual arousal by two variables, $u(t)$ and $v(t)$, which are taken to represent current levels of physiological and psychological/cognitive arousal, respectively. In doing so, we address some of the above-mentioned criticisms of the Masters–Johnson model in terms of not appropriately accounting for the psychological component of the arousal.

^{14–16}The equation for the physiological component of arousal $u$ is taken to be of the form

With the average duration of a sexual intercourse being estimated to be around 5.4 min,^{36} which includes between 100 and 500 thrusts,^{37} instead of considering physical stimulation as a periodic forcing $ E ucos\u2061(\omega t)$ with some very high frequency $\omega $, we rather represent it as a constant $ E u$ over the entire duration of the sex cycle.

The term $(\u2212v)$ on the right-hand side of Eq. (6) may seem counter-intuitive, because it appears to suggest that the growth of physiological arousal decreases with increasing psychological arousal, and hence, it requires some justification. The first reason for choosing this term comes from PET and fMRI studies of men, which showed that an orgasm was associated with a decreased activity of amygdala^{38} and temporal lobe,^{39} as well as with deactivation of frontal cortical regions^{22,40,41} and orbitofrontal cortex.^{22} This suggests that reaching an orgasm can be characterized as psychologically “letting go,”^{42} and having a term that decreases with psychological arousal achieves exactly that. The second reason is that since a psychological response to sex has been shown to be neurologically similar to response to other pleasures,^{21,43,44} one can rely on known results suggesting that optimal performance is associated with some intermediate level of mental/psychological arousal. An inverted-U curve of performance and arousal describing this is often associated with the so-called Yerkes–Dodson law in the context of reward-based training,^{45} which has been subsequently applied in a variety of other contexts looking at dependence of performance on arousal.^{46–48} A simple intuitive explanation behind this observation is that low levels of mental arousal are associated with boredom/apathy, whereas very high levels result in anxiety, and hence, it is the intermediate levels of psychological arousal that result in the optimal level of physical performance. In Fig. 3(b), we have plotted the nullcline associated with the right-hand side of Eq. (6), where the contribution of psychological arousal is given either by the term $(\u2212v)$ as stated in that equation or, alternatively, as a cubic function of the form $v\u2212 v 3/3$, which could represent the growth of physiological arousal for small values of psychological arousal, and a reduction in physiological arousal for high values of psychological arousal. This cubic function represents the Yerkes–Dodson law in terms of there being an optimal level of physiological arousal that provides the highest level of (physiological) performance. We note from this figure that while the specific values of variables on the nullclines have changed, its qualitative shape is the same as that for the linearly decreasing function of the form $(\u2212v)$. This means that the structure of steady states and their stability would not change if we were to consider the more complicated cubic function. Hence, to simplify the analysis, we will use the functional form given in Eq. (6).

^{21,43,44}as well as from the so-called

*opponent process theory*.

^{49–51}This psychological and neurological theory suggests that the emotional process can be split into two parts: the A-process, which is a fast response to an external stimulus [the first two terms in Eq. (8)], and the B-process that is responsible for returning to the state of psychological homeostasis. This theory has already been effectively used to study a number of psychological processes, such as drug addiction and withdrawal,

^{52–55}as well as sleep.

^{56,57}Finally, whereas physiological responses are rather quick and are controlled by the autonomic nervous system, processing of psychological stimuli is (much) slower; hence, we introduce a scaling parameter $0<\u03f5\u226a1$. In both Eqs. (6) and (8), parameters $ E v 0$, $a$, and $b$ are assumed to be positive, while external stimulations $ E u$ and $ E v$ are non-negative.

## IV. STABILITY ANALYSIS AND NUMERICAL SIMULATIONS

When $ u \u2217$ is lying on the middle branch, i.e., $ u A< u \u2217< u B$, the condition of uniqueness of the steady state $ E \u2217$ guarantees that $det( J 0)=\u03f5 [ a \u2212 b f \u2032 ( u \u2217 ) ]>0$. Then, the steady state $ E \u2217$ is stable for $ f \u2032( u \u2217)\u2212b\u03f5<0$, unstable for $ f \u2032( u \u2217)\u2212b\u03f5>0$, and undergoes a supercritical Hopf bifurcation when $ f \u2032( u \u2217)\u2212b\u03f5=0$. We note that since $a\u2212b f \u2032(u)>0$ for all values of $u$, this means that $det( J 0)=\u03f5 [ a \u2212 b f \u2032 ( u \u2217 ) ]$ is always positive, and hence, $\lambda =0$ is not a root of the characteristic equation. This then implies that the only change of stability can occur through a Hopf bifurcation, where a pair of complex conjugate eigenvalues cross the imaginary axis. Similarly to the case mentioned earlier for the left branch, at $ [ f \u2032 ( u \u2217 ) + b \u03f5 ] 2\u22124\u03f5a=0$, the two unstable eigenvalues will coincide on the real axis, and an unstable focus will become an unstable node. As we move further along the middle branch to the right, at some point, the unstable node will turn back into an unstable focus, and eventually, the inverse supercritical Hopf bifurcation will take place, where the steady state $ E \u2217$ will regain its stability. Finally, on the right branch $ u \u2217> u B$, the steady state $ E \u2217$ is stable. The points $ u \u2217= u A$ and $ u \u2217= u B$ are always stable, and they are either stable nodes if $ b 2\u03f5\u22654a$ or stable foci otherwise. In light of condition (10) that ensures uniqueness of the steady state $ E \u2217$ and the fact that $0<\u03f5\u226a1$, in most cases, the points $ u \u2217= u A$ and $ u \u2217= u B$ will be stable foci.

In Fig. 5, we present one-dimensional bifurcation diagrams obtained by fixing all of the parameters except one, and then for each value of bifurcation parameter, we note the corresponding steady-state value of psychological arousal $ v \u2217$, as well as stability type of the steady state $ E \u2217$. Figure 5(a) illustrates the sequence of transitions shown in Fig. 4 when the level of physiological stimulation $ E u$ is fixed, while the level of psychological stimulation $ E v$ is increasing. This corresponds to moving $v$-nullcline upward, while $u$-nullcline remains the same. In this case, we observe transitions from an unstable node to an unstable focus, followed by an inverse supercritical Hopf bifurcation, resulting in the stable focus followed by the stable node. Increasing the level of physiological stimulation $ E u$ does not alter the shape of the bifurcation curve but results in higher values of $ v \u2217$ and stability changes occurring at slightly higher values of $ E v$. Figure 5(b) illustrates a similar sequence of transitions, though starting all the way on the right branch of the $u$-nullcline and moving to the left along this nullcline for increasing values of $a$.

Figure 6(a) illustrates the phase plane of the system (6)–(8). Similarly to the well-known case of a FitzHugh–Nagumo model, the left and right branches of the $u$-nullcline are attracting, while the middle branch is repelling, and there exists a separatrix shown in green such that initial conditions that lie in the phase space to the left of this separatrix return to the steady state $ E \u2217$ without ever reaching the right branch of the $u$-nullcline. We define $u= u B$ as the critical point in the sense that if the level of physiological arousal has exceeded this value before returning to $ E \u2217$, this will be interpreted as an orgasm having happened. Numerically, separatrix can be found by either taking initial conditions on the right branch of $u$-nullcline and integrating the system backward in time^{58} or by taking the single point $(u,v)=( u B,f( u B)+ E u)$ as an initial condition and integrating the system backward in time. Importantly, in the same way as with the FitzHugh–Nagumo model, there is no clearly defined threshold in the model, and hence, separatrix plays the role of a quasi-threshold manifold.^{59} Thus, separatrix has an important biological interpretation as a manifold separating orgasmic and anorgasmic trajectories, which are illustrated in the same plot in red and blue. These trajectories are reminiscent of “canard” trajectories that are found in a number of slow–fast systems, including the FitzHugh–Nagumo model.^{60,61} Canards are characterized by following a slow manifold for a long period of time before making a fast excursion in the phase space. Such solutions often appear as a result of folded nodes and singular Hopf bifurcations,^{60,61} and they owe their name to the French school of mathematicians,^{62–64} who noted a similarity between the shape of canard cycles in the phase plane and ducks (“canard” is a duck in French). In our case, orgasmic trajectory can be characterized as a “canard with head,” while an anorgasmic trajectory is a “canard without head.”^{65,66}

Temporal dynamics of orgasmic (red) and anorgasmic (blue) trajectories is shown in Figs. 6(b) and 6(c). They reproduce the structure of the Masters–Johnson sexual response cycle by going sequentially through EPOR stages. The cycle starts with a fast initial excitation, followed by a steady, but slower growth of arousal during the plateau stage, and then, depending on the initial conditions and the values of parameters, the trajectory will either cross the orgasmic threshold, and then an orgasm will be followed by a resolution, or resolution will occur without an orgasm being achieved. Unlike the original Masters–Johnson theory, which assumed arousal to stay constant during the plateau phase, in our model, it continues to grow, albeit at a slower rate, which agrees with several studies.^{17,67} Looking more carefully at an orgasmic trajectory, after reaching an orgasm, it will return to a steady state $ E \u2217$, but since this steady state is a stable focus, approach toward it will occur in an oscillatory manner. This agrees with an observation that an orgasm (both in males and females) is characterized by regular periodic pelvic muscle contractions.^{7,68,69} Furthermore, if we look more carefully at the orgasmic trajectory shown in red, we notice that straight after an orgasm, for a short period of time, this trajectory has negative values of physiological arousal. This fits with a clinical observation that due to penile hyper-sensitivity after an orgasm, ongoing physical stimulation is often felt as unpleasant or even painful,^{70,71} extending for some part of the refractory period. Physiologically, it has been suggested that this happens due to a substantial release of prolactin immediately after an orgasm, which produces the feeling of sexual satiety, while also causing a notable drop in the dopamine level resulting in feeling tired.^{72,73}

## V. STOCHASTIC EFFECTS

Since people are exposed on a daily basis to a variety of physical and emotional stimuli, not to mention circadian and other changes of hormonal levels, when studying the dynamics of such complex processes as a sexual response, it is important to consider possible effects of stochasticity. Data show that men experience an average of 11 erections during the day, plus 3–5 instances of nocturnal penile tumescence, usually associated with a period of REM-sleep.^{75} Since majority of these are not caused by any organized sexual activity, they can be considered stochastic perturbations of an otherwise (deterministically) stable state. Understanding statistical properties of spontaneous erections is also important from the perspective of using monitoring of nocturnal penile tumescence and rigidity (NPTR) in clinical practice for a differential diagnosis of a psychogenic vs physiological (organic) erectile dysfunction.^{76–78}

^{79,80}

### A. Stochastic oscillations near equilibrium

^{81–83}In this case, the covariance matrix $\Sigma $ for fluctuations is constant and is given by the solution of the Lyapunov equation,

^{83}known in the context of control problems as the continuous-time algebraic Riccati equation,

^{84}

^{,}

^{81,83}

^{81,83}

^{85–87}Using SSF, we can find the so-called

*confidence ellipse*

^{88}

^{,}

*stochastic amplification*,

^{89,90}where on average, the system exhibits decaying oscillations toward deterministically stable steady state $ E \u2217$, whereas in individual stochastic realizations, it still shows sustained stochastic oscillations. Physiologically, this is a really important observation, as it highlights the fact that while generally being at rest, small physiological and/or psychological excitations can indeed result in periods of spontaneous penile tumescence, as confirmed by clinical observations.

^{75}

*coherence*that can be defined as follows. We consider an interval of frequencies $[ \omega 1, \omega 2]$ around the dominant frequency $ \omega 0$, i.e., $ \omega 1< \omega 0< \omega 2$, and compute a quantity (cf. expressions in Eqs. (2.8) and (2.9) in Ref. 90)

^{90}

Figure 8 illustrates spectra of stochastic oscillations, as given by (23), together with their variance and coherence. We observe for fixed values of all other parameters, increasing the strength $ E v$ of psychological stimulation results in making the frequency distribution of stochastic oscillations much closer to uniform, effectively smearing out the dominant frequency. In contrast, increasing the rate $b$ of relaxation of psychological arousal to its state of homeostasis makes the frequency distribution sharper around the dominant frequency. Reducing psychological stimulation $ E v$ or reducing the rate $a$, at which psychological arousal increases with physiological stimulation, leads to a higher variance of stochastic oscillations around the deterministically stable steady state $ E \u2217$, as well as a higher level of coherence of these oscillations, reaching its maximum value of 1 at the boundary of Hopf bifurcations, where this steady state loses its stability, giving rise to deterministically stable periodic solutions.

### B. Large deviations: An orgasm as an escape

*Large deviation theory*(LDT)

^{91}provides a quantitative characterization of this decay of probabilities. If we introduce an action functional, also known as the rate function, as

*quasi-potential*$V( x 1, x 2)$

^{92–94}that characterizes long-term behavior of the system as

^{91}

^{,}

*instanton*,

^{79,92,95}is effectively the least unlikely trajectory that escapes $ B( E \u2217)$, and it produces the smallest action $ S T[ x \u2217(t)]$; i.e.,

^{91}

^{,}

^{79,92}

^{59}MPEP can be viewed as solutions of the extended system that originate in the small neighborhood of the saddle $( E ~ \u2217, 0)$ and cross the separatrix, indicating escape toward the right branch of the $u$-nullcline.

One method of computing MPEP is the method of an action plot,^{96} in which a large number of initial conditions are chosen to be uniformly distributed on a small circle around the steady state $( E ~ \u2217, 0)$, and then for each of them, the Hamiltonian system (33) is integrated until the solution crosses the separatrix. At that point, the value of action $S$ is recorded, and the trajectory associated with the minimum value of the action is taken to represent an optimal escape path. Another approach is the geometric minimum action method (gMAM),^{92–94} in which the fact that the action $S( x)$ is independent of path parameterization allows one to replace that action by an action that is parameterized by normalized arc length instead of time. This allows one to overcome the numerical difficulty associated with the fact that since the escape trajectory starts at the steady state, time $T$ in the definition of action is infinite, whereas it would be replaced by a finite arc length once action is replaced by a geometric action. Several versions of the gMAM method have been proposed, including the simplified gMAM^{92} and the adaptive minimum action method,^{97,98} and they have been successfully used to compute MPEP for a number of systems, including the FitzHugh–Nagumo model,^{59} a system with overdamped double-well potential,^{99} a Ludwig spruce budworm model,^{79} a predator–prey model with migration,^{92} and a Maier–Stein model^{100} that has long served as a benchmark for LDT calculations. Below, we will use gMAM to compute the optimal path for escape from the saddle $( E ~ \u2217, 0)$ to points on a separatrix.

^{101}In this respect, a more comprehensive description is provided by the global quasi-potential as computed over an entire region of interest. A first-order accurate ordered upwind method (OUM)

^{102}and an ordered line integral method (OLIM)

^{101,103}have been proposed for numerically computing quasi-potential as a solution of the Hamilton–Jacobi equation associated with the Hamiltonian (32). OUM has been used to find quasi-potentials and to compute MPEP in the Higgins model for glycolysis,

^{104}the Zeldovich–Semenov model for an exothermic reaction in a continuously stirred tank reactor,

^{105}and the FitzHugh–Nagumo model of neural excitation,

^{59}among others. Once the quasi-potential $ U E \u2217$ is found, MPEP that goes from the attractor $ E \u2217$ to a given point $ x$ can be found by integrating the equation

^{101–103}

To compute quasi-potential and MPEP for the system (26), we considered a region $[\u22121.5,7]\xd7[0,11]$ that was divided by a rectangular mesh $2000\xd72000$, over which quasi-potential was computed in MATLAB with a modification of the OUM code developed earlier for the Higgins model.^{106} The resulting quasi-potential and its contour lines over the $ u ~$– $ v ~$ plane are shown in Figs. 9(a) and 9(b). The shape of quasi-potential can be described as a deep well around the steady state, with a plateau-like area further away. Since we are interested in the MPEP, where stochastic trajectories escape a stable steady state by crossing a quasi-threshold associated with the separatrix, we investigated how the value of quasi-potential changes along this separatrix. This analysis reveals that the change is non-monotonic, and there is a point that we denote as $ S 0$ that represents the minimum of quasi-potential along the separatrix. Since, according to (30), such a point would provide the shortest escape time and would also correspond to the optimal escape path,^{91} we choose this point as an initial condition for shooting backward to the steady state using the flow along the quasi-potential (34), as well as for gMAM calculation.

For an action plot, we took 20 000 initial conditions chosen as $ x 0 i= E \u2217+\delta x 0 i$, where $\delta x 0 i= ( R cos \u2061 \varphi i , R sin \u2061 \varphi i ) T$, $R=0.1$, $i=1\u202620,000$, uniformly distributed around the circle of radius $R$ centered on the steady state, with initial conditions for conjugate momenta $ p$ and action $S$ chosen accordingly.^{96,107} Then, we solved the Hamiltonian system (33) for each of these initial conditions until the solution reached the separatrix, at which point the value of action was recorded. The resulting action plot is illustrated in Fig. 9(c), and it has the shape similar to that for other models;^{96,104} namely, there is a minimum value of action, and there are several other initial conditions that yield values of action very close to the minimum. The value of action at the point $ S 0$, associated with the minimum of quasi-potential along the separatrix, differs from the overall minimum value of action by less than $0.04%$, and hence, this point $ S 0$ can be chosen as the target point for computing MPEP in terms of minimizing action, as computed using an action plot. Hence, we chose the phase $\varphi =0.0014646947$ to set up initial conditions and then computed MPEP by following the Hamiltonian system (33), as shown in Fig. 9. Having $ E \u2217$ and $ S 0$ as two end points, we employed gMAM to compute MPEP, which is illustrated in blue in the same plot, and it closely follows the trajectory computed using the action plot. Finally, we also show in this figure the optimal path as described by (34), where we used quasi-potential that had already been computed with OUM.

In terms of interpretation of these results for our model, we note that the point $ S 0$ that corresponds to the minimum of quasi-potential along the separatrix is characterized by the values of physiological and psychological arousal that do not significantly exceed their values at the steady state, with a larger difference being observed for physiological arousal. Importantly, these results not only show that it is possible to achieve an orgasm as a result of stimuli that arise from small amounts of stochastic physiological and/or psychological stimulation, but they also suggest that actually when this happens, with it being a stochastic escape through a separatrix, the levels of physiological and psychological arousal do not need to be high. One relevant example from clinical practice concerns males with spinal cord injuries (SCIs), who are known to experience difficulties with achieving an orgasm post-injury.^{108–110} While the reasons for this can be both psychogenic and organic, it has been shown that with damage to the spinal cord, there may be greater reliance on psychogenic stimulation for achieving arousal and orgasm,^{111,112} and while the distinction between reflexive and psychogenic erection in able-bodied males may be less clear,^{113} for males with SCI, this difference is crucial, and the possibility of stochastic escape to an orgasm even in the absence of a sufficiently high level of physiological arousal may prove very important.

The idea of small-noise stochastic escape from the basin of attraction of a stable steady state is that the probability of escape is very low and decreasing with noise magnitude, which, in practical terms, means that stochastic escape is indeed possible, though this happens very rarely. A good proxy for stochastic escape associated with an orgasm can be provided by spontaneous nocturnal emissions, which are ejaculations that happen during sleep as an autonomous reflex mediated by sympathetic nervous systems without any stimulation.^{114,115} While orgasm and ejaculation are two separate physiological processes and it is possible for males to experience an orgasm in their sleep without ejaculation,^{5} though in the majority of cases, both of these processes happen at the same time.^{75} Due to the sensitive nature of the subject, it is difficult to obtain robust estimates of the frequency of nocturnal emissions, with significant variations observed for different age groups (men aged 19–20 having the highest frequency^{5}) and different societal groups.^{116} A more recent survey^{114} estimated the frequency to range from once a month for males in their late teens and early 20s, to once or twice per year for older adult males. In almost all cases, males reported having a sex dream when they experienced a nocturnal emission,^{117} which can be interpreted within the framework of our model as stochastic escape from the basin of attraction of an otherwise stable baseline state of low levels of physiological and psychological arousal under the influence of a small amount of stochastic psychological stimulation.

## VI. DISCUSSION

In this paper, we have studied a mathematical model of a sexual response in human males with account for physiological and psychological arousal. Depending on the values of system parameters and initial conditions, the model is able to exhibit solutions, whose dynamics reproduces that of the Masters–Johnson sexual response cycle, and more specifically, progression through the stages of excitation–plateau–orgasm–resolution, or the same type of progression without achieving an orgasm. In the absence of a fixed threshold in the model, dynamics in the phase space is mediated by the presence of a separatrix, which plays the role of a quasi-threshold, with only solutions that cross the separatrix proceeding to reach the right branch of the corresponding nullcline, which is associated with achieving an orgasm. Dynamically, such trajectories have the behavior similar to that of canard trajectories in other slow–fast system in the sense that they follow an unstable slow manifold for an extended period of time, before making a fast excursion in the phase space ending in a stable steady state. Numerical solutions illustrating orgasmic trajectories also exhibit a short period of negative physiological arousal immediately following an orgasm, a phenomenon well known from clinical practice.

To obtain a better understanding of the effect stochastic perturbations may have on system dynamics, we have analyzed stochastic fluctuations in the neighborhood of a stable steady state under the influence of noise in both types of arousal, while also ensuring that noise terms are acting on their respective time scales. The results indicate the possibility of stochastic amplification, where sustained stochastic oscillations around a deterministically stable steady state are observed, and we have computed their spectrum, variance, and coherence. We have also looked at confidence domains that delineate regions in the phase space surrounding that steady state, where stochastic solutions are contained with a certain probability. Since an orgasm is associated with a significant increase in physiological arousal compared to baseline, we have employed large deviation theory to compute quasi-potential, as well as optimal escape paths from the neighborhood of the stable steady state associated with baseline arousal. These calculations show that while rare, under the influence of small stochastic perturbations, it is possible for trajectories to escape basin of attraction of this stable steady state, cross the separatrix, and proceed for a large excursion in the phase space associated with achieving an orgasm, which again is a clinically known phenomenon of nocturnal emissions often caused by dreams of sexual nature.

Despite its relative simplicity, the model proposed and studied in this paper is able to qualitatively reproduce all essential features of a sexual response in human males, as observed in clinical practice. There are a number of directions, in which the work presented in this paper could be further extended. One interesting problem would focus on more detailed analysis of dynamics in the neighborhood of a baseline state. Masters and Johnson^{7} were the first to note that pelvic muscle contractions during an orgasm appeared to have the period of around 0.8 s almost universally in males and females. Since our results provide analytic expressions for the frequency of decaying oscillations toward this state, as well as their decay rate, comparing these against clinical measurements could provide estimates of some of the fundamental parameters characterizing physiology of a sexual response. Additionally, spectra of stochastic oscillations can be compared to measurements of NPTR, which could provide some understanding of the structure of NPTR as observed in people with normal sexual function and those with erectile dysfunction. This could then help produce prognostic criteria for a better diagnosis of psychogenic vs physiological erectile dysfunction, thus improving clinical monitoring and treatment of this condition, especially with recent advances in wearables and medical monitoring devices.

While this paper has been concerned with the dynamics of aggregated levels of physiological and psychological arousal, another very interesting potential avenue of research concerns a more detailed representation of neural dynamics associated with a sexual response. An orgasm has been called “an altered state of consciousness,”^{118} and several neurological studies have shown that neurons in some parts of the brain exhibit the level of synchrony during an orgasm, which is very similar to that observed during epileptic seizures.^{119–121} The significance of this observation is that one can gain insights into emergence and control of neural synchrony during a sexual response and an orgasm by using mathematical techniques developed for analysis of various phenomena in systems of coupled neurons, such as amplitude/oscillation death,^{122} chaotic synchronization,^{123} and chimera states.^{124} In this context, it has been suggested that periodic physiological stimulation, which we modeled as a constant term due to its high-frequency nature when considered on the timescale of a sexual response cycle, can actually result in the entrainment of oscillations among groups of neurons,^{119} thus suggesting the possibility of another feedback mechanism between physiological and psychological components of the arousal.

One could also revisit the relation between arousal and performance as used in the equation for physiological arousal and reconsider it from the perspective of reversal theory.^{125,126} In this framework, depending on the tone of activity, when in the playful state, the increase in psychological arousal moves one from boredom to excitement, whereas in the serious state, a similar increase in psychological arousal moves between relaxation and anxiety. This allows a more nuanced representation of interactions between the level of psychological arousal and physiological pleasure, as is associated with a sexual response, not to mention that it can provide a wider range of scenarios for switching between playful and serious states, and the corresponding changes in arousal. When combined with the above-mentioned data on neural correlates of a sexual response and an orgasm, this could yield a much more accurate representation of a human sexual response. Finally, with insights obtained from the results in this paper, the next major step would be to consider how one could modify modeling framework to study the sexual response in females, which is known to be much more complex from the perspective of interactions between physiological and psychological arousal, as well as in terms of a much wider range of different temporal dynamics that it can exhibit.^{7,8} More specifically, one could incorporate in the model additional feedback mechanisms proposed by Basson^{8} in her circular model of a female sexual response, which pointed to desire and excitation arising as a by-product of stimulation in addition to spontaneous desire associated with stochastic inputs, and the role of sexual motivation as a factor influencing the degree of receptiveness to stimulation.

## SUPPLEMENTARY MATERIAL

See the supplementary material for bifurcation analysis and numerical simulations of the model with a physiological response to psychological arousal being represented by a cubic function modeling the Yerkes–Dodson law.

## ACKNOWLEDGMENTS

The authors would like to thank Dr. Yang Li for help with numerical codes used for calculation of quasi-potential. The authors acknowledge financial support from the Dr Perry James (Jim) Browne Research Centre.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**K. B. Blyuss:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Visualization (equal); Writing – original draft (equal). **Y. N. Kyrychko:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Visualization (equal); Writing – original draft (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

## REFERENCES

*Mathematical Physiology II: Systems Physiology*

*Applied Mathematical Models in Human Physiology*

*An Introduction to Mathematical Physiology*, University of Oxford Lecture Notes on Mathematical Physiology (University of Oxford, 2021); see https://courses.maths.ox.ac.uk/pluginfile.php/12781/mod_resource/content/1/printed_notes/physiolbook.pdf.

*Sexual Behavior in the Human Male*

*Sexual Behavior in the Human Female*

*Human Sexual Response*

*Appetite-Neural and Behavioural Basis*, edited by C. R. Legg and D. A. Booth (Oxford University Press, Oxford, 1994), pp. 127–164.

*Human Sexuality and Its Problems*

*Handbook of Sexology*, edited by J. Money and H. Musaph (Elsevier, Amsterdam, 1977), pp. 327–350.

*Alternative Approaches to the Study of Sexual Behavior*, edited by D. Byrne and H. Kelly (Erlbaum, Hillsdale, NJ, 1986), pp. 1–12.

*The Textbook of Clinical Sexual Medicine*, edited by W. W. IsHak (Springer, 2017), pp. 39-50.

*3rd International Symposium on Systems Research, Informatics and Cybernetics*(Baden-Baden, Germany, 1991).

*The Proceedings of the 2000 International Conference on Mathematics and Engineering Techniques in Medicine and Biological Sciences*(METMBS2000) (Las Vegas, NV, 2000).

*The Use of Models in the Social Sciences*, edited by L. Collins (Routledge, New York, 1975).

*Catastrophe Theory—Selected Papers 1972-1977*

*Nonlinear Ordinary Differential Equations—An Introduction for Scientists and Engineers*

*Patterns of Sexual Arousal: Psychophysiological Processes and Clinical Applications*

*et al.,*

*User Modeling, Adaptation, and Personalization*, edited by S. Carberry, S. Weibelzahl, A. Micarelli, and G. Semeraro (Springer, Rome, 2013), pp. 127–138.

*Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting*

*et al.,*

*Handbook of Sexual Dysfunction*, edited by R. Balon and R. T. Segraves (CRC Press, Boca Raton, FL, 2005).

*The Orgasm Answer Guide*

*Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences*

*Algebraic Riccati Equations*

*Random Perturbations of Dynamical Systems*

*Recent Progress and Modern Challenges in Applied Mathematics, Modeling and Computational Science*, Fields Institute Communications, edited by R. Melnik, R. Makarov, and J. Belair (Springer, New York, 2017), pp. 17–55.

*The Psychobiology of Consciousness*, edited by J. M. Davidson and R. J. Davidson (Plenum Press, New York, 1980), pp. 271–332.

*Synchronization: A Universal Concept in Nonlinear Sciences*

*The Dangerous Edge: The Psychology of Excitement*

*Zigzag: Reversal and Paradox in Human Personality*