Sex, ducks, and rock "n" roll: mathematical model of sexual response

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.

Mathematical models have proven to be very effective tools for understanding the dynamics of various physiological processes and have significantly contributed to improving monitoring and treatment for a number of diseases. Of course, success of these models strongly relies on the possibility of direct experimental measurements, as well as availability of reasonably accurate experimental data that are used for model derivation. One exception is a physiological description of human sexual function, which has been largely based on self-reporting surveys, and until very recently not amenable to accurate measurements in controlled conditions, which prevented mathematical analysis of this process. In this work, we propose a mathematical model of a sexual response cycle in the human male, which is based on the ideas of the classical Masters-Johnson theory with a number of modifications reflecting subsequent improvements of this theory. Formulating the model as an excitable slow-fast system of two equations for the levels of physiological and psychological arousal, we explore the roles played by excitations applied to either of these components and perform numerical simulations to illustrate how the model can exhibit different types of dynamics observed in clinical practice. To account for the fact that people can receive a variety of stimuli, we analyze the effects of stochasticity on model dynamics by studying the structure of stochastic oscillations close to an equilibrium. Using large deviation theory, we are able to find optimal stochastic escape paths that show how a sexual response progresses toward an orgasm under the influence of small stochastic perturbations.

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][2][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 → plateau → orgasm → 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][9][10][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][15][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][21][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
Following the development of catastrophe theory by Thom and Zeeman in the 1960s and 1970, it has subsequently been applied in a variety of contexts, ranging from ecology 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ẋ with Steady states of this model satisfy a cubic equation and depending on the values of parameters a and b, there can be one or three possible roots of this equation. Figure 1(a) shows the surface of steady states in the space of parameters, which illustrates the 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. Motivated by ideas from the studies of periodically forced systems, Hubey 29,30 suggested that a sexual response can be mathematically modeled as a Duffing oscillator, where θ is some characterization of sexual arousal, and the parameters are described in terms of mechanical interpretation: α is a damping coefficient, β controls linear stiffness, γ represents nonlinearity of restoring force, and f and are, respectively, amplitude and frequency of forcing. Expression for this amplitude response curve can actually be found analytically, using the method of harmonic balance, which gives 34 This amplitude response exhibits hysteresis: if one gradually increases forcing frequency , this results in the initial increase in the amplitude of oscillations followed by a drop onto a lower branch, where a further increase in the forcing frequency would lead to a decrease in the amplitude of oscillations. Going in the opposite direction, starting with a high forcing frequency and decreasing it, there would be a jump onto the upper branch of the amplitude response curve. Hubey 29,30 plotted an amplitude response curve for the model (2), i.e., amplitude of oscillations |θ | around zero steady state of the unforced model, as a function of forcing frequency , 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 twodimensional 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
Following an observation 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][15][16] The equation for the physiological component of arousal u is taken to be of the formu where f(u) is an inverse N-shaped function, as illustrated in Fig. 3(a), the term (−v) describes increasing physiological arousal by means of reducing the level of psychological arousal, and E u is the level of applied physical stimulation. Function f(·) is chosen to be of the following form that is inspired by earlier works on excitable systems: with two critical points located at u A and u B and sufficiently smooth functions f L (u), f M (u), and f R (u) describing left, middle, and right branches, respectively. In general, one can choose three different functional forms for f L (u), f M (u), and f R (u), but as long as they satisfy matching conditions = 0, they will still provide a continuously differentiable right-hand side of Eq. (6). We will choose these functions to be In the particular case where f(u) = u − u 3 /3, we recover the classical FitzHugh-Nagumo-like function. The purpose of choosing a higher degree of u on the right branch of function f(u) is to provide a mechanism for a faster resolution phase after an orgasm. 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 u cos(ωt) with some very high frequency ω, we rather represent it as a constant E u over the entire duration of the sex cycle.
The term (−v) 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][47][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 (−v) as stated in that equation or, alternatively, as a cubic function of the form v − 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 (−v). 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).

Dynamics of a psychological component of arousal is modeled by the following equation
where E v0 is a baseline level of psychological arousal, E v is external psychological stimulation associated with visual, verbal/auditory, and other stimuli, and it is assumed that as physiological arousal grows, this will, in turn, also increase psychological arousal, as represented by the term au. Justification for the last term (−bv) comes from earlier studies that have demonstrated a similarity between a psychological response to sex and that to other pleasures, 21,43,44 as well as from the so-called opponent process theory. [49][50][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 Bprocess 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 < 1. In both Eqs. (6) and (8), parameters E v0 , a, and b are assumed to be positive, while external stimulations E u and E v are non-negative.

ARTICLE scitation.org/journal/cha
These nullclines are illustrated in Fig. 3(a), and we observe that increasing/decreasing either of E u or E v results in moving the corresponding nullcline upward or downward, which then affects, where the two nullclines intersect. We will assume that the following condition holds: which guarantees that the model (6)-(8) always has a single steady state E * = (u * , v * ). Introducing auxiliary parameters k = a/b and d = (E v0 − E v ) /a, we then have the following classification: E * lies on the right branch.
For the particular choice of the cubic function f given in (7), the locations of its two critical points are given by Furthermore, since this function is cubic, it will have an inflection point where f (u infl ) = 0 and f (u infl ) reaches its maximum; therefore, the condition (10) for uniqueness of the steady state turns into Jacobian of linearization near the steady state E * is given by with Tr(J 0 ) = f (u * ) − b and det(J 0 ) = a − bf (u * ) . If the steady state E * lies on the left branch, i.e., u * < u A , we have f (u * ) < 0, and, therefore, for positive values of parameters, Tr(J 0 ) < 0 and det(J 0 ) > 0, which implies that the steady state When u * is lying on the middle branch, i.e., u A < u * < u B , the condition of uniqueness of the steady state E * guarantees that det and undergoes a supercritical Hopf bifurcation when f (u * ) − b = 0. We note that since a − bf (u) > 0 for all values of u, this means that det(J 0 ) = a − bf (u * ) is always positive, and hence, λ = 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 (u * ) + b 2 − 4 a = 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 * will regain its stability. Finally, on the right branch u * > u B , the steady state E * is stable. The points u * = u A and u * = u B are always stable, and they are either stable nodes if b 2 ≥ 4a or stable foci otherwise. In light of condition (10) that ensures uniqueness of the steady state E * and the fact that 0 < 1, in most cases, the points u * = u A and u * = u B will be stable foci.
In order to perform numerical bifurcation analysis and simulations, we fix the coefficients in the definition of function f(u) to be and the corresponding bifurcation diagram for the steady state E * is shown in Fig. 4. We observe that for a fixed value of physiological stimulation E u , increasing the level of psychological stimulation E v results in the steady state E * changing its stability from an unstable node to an unstable focus by means of the above-mentioned collision of characteristic eigenvalues on the real axis, followed by a transition to a stable focus via an inverse supercritical Hopf bifurcation, and then to a stable node, after characteristic eigenvalues collide again on the real axis. If one fixes the value of psychological stimulation and increases physiological stimulation, the steady state E * experiences the same sequence of stability changes, but in the reverse order. When looking at stability of E * in terms of parameters a and b, as shown in Fig. 4(b), we observe that there is a bound on the values of parameter b as given by the condition a − bf (u) > 0, which ensures the existence of the steady state E * . As the value of parameter a that quantifies how psychological arousal grows with increasing physiological arousal, we observe the same sequence of transitions, as if we were increasing physiological stimulation E u . For sufficiently high values of a, most of the a-b parameter plane corresponds to the steady state E * being a stable focus. However, there is an intermediate range of a values, for which, provided the value of parameter b, which characterizes the rate of return to homeostasis in the level of psychological arousal, is sufficiently small, the steady state E * is a stable node. 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 steadystate value of psychological arousal v * , as well as stability type of Colors denote a stable node (black), a stable focus (blue), an unstable focus (red), and an unstable node (yellow).
the steady state E * . 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 * and stability changes occurring at slightly higher values of E v .   (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 * 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 * , 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 Chaos ARTICLE scitation.org/journal/cha 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 * , 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. 74 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. [75][76][77]

ARTICLE scitation.org/journal/cha
In light of these observations, we modify the model (6)-(8) by writing it as a system of Langevin equations, where D 1 and D 2 are intensities of noise in physiological and psychological components of arousal, respectively, and η 1 and η 2 are uncorrelated Gaussian white noises with η i = 0 and η i (t)η j (t ) = δ ij δ(t − t ) for i, j = 1, 2. An extra factor √ in the noise term in the second equation of stochastic model (13) has been introduced to ensure noises act on their respective timescales. 78,79 A. Stochastic oscillations near equilibrium As mentioned earlier, for most of the time, the level of sexual arousal, both physiological and psychological, is maintained at some constant small level, which, in the deterministic model, corresponds to a stable steady state E * . To explore the emergence and statistical properties of arousal in the neighborhood of E * , we linearize the model (13) near this deterministically stable steady state and obtain a system of equations for stochastic fluctuations ξ 1,2 (t), where J 0 is the Jacobian (12) evaluated at E * , and we have the following connection to original variables: u(t) = u * + ξ 1 (t) and v(t) = v * + ξ 2 (t). Since we are considering fluctuations associated with linearization near a stable steady state E * , which is characterized by a pair of complex conjugate eigenvalues with negative real parts, the system (14) describes a two-dimensional Ornstein-Uhlenbeck process. [80][81][82] In this case, the covariance matrix for fluctuations is constant and is given by the solution of the Lyapunov equation, 82 known in the context of control problems as the continuous-time algebraic Riccati equation, 83 where D is the diffusion matrix, The solution of this equation can be found in a closed form, 80,82 This gives the variance of stochastic fluctuations in the level of physiological arousal as The Fokker-Planck equation for the probability density ρ(ξ , t) associated with system (14) has the form Due to the fact that eigenvalues of the Jacobian J 0 have negative real parts, this equation has a solution in the form of stationary probability density, 80,82 Since the covariance matrix describes the magnitude and the spatial arrangement of stochastic fluctuations around the exponentially stable steady state E * , for this reason, it is also known as the stochastic sensitivity function (SSF). [84][85][86] Using SSF, we can find the so-called confidence ellipse where k 2 = − ln(1 − p) and p is a fiducial probability; i.e., random solutions lie in the interior of this ellipse with probability p. Using eigenvalues (µ 1 , µ 2 ) and corresponding eigenvectors (v 1 , v 2 ) of the covariance matrix , the confidence ellipse can be expressed in the form 87 with new variables ζ 1,2 = (x − x 0 ) T v 1,2 . Figure 7(a) illustrates the confidence ellipse in the u-v plane around the steady state E * for different levels of noise. In Fig. 7(b), we show a comparison of the deterministic trajectory with a single stochastic realization of the system (13), illustrating the phenomenon of stochastic amplification, 88,89 where on average, the system exhibits decaying oscillations toward deterministically stable steady state E * , 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. 74 While variance provides some quantitative characterization of stochastic fluctuations, further insights into their structure can be obtained by computing their spectra. To this end, we take a Fourier transform of the system (14), which yields This then gives the power spectrum matrix as where † denotes a Hermitian conjugate, and, for notational convenience, we have denoted as a ij coefficients of the Jacobian matrix J 0 = a 11 a 12 a 21 a 22 .
For stochastic fluctuations in physiological arousal, we then have with where all four parameters A, B, C, and D are positive. Differentiating S u as a function of z, we find that S u has extrema at and since we require z > 0 due to the relation z = ω 2 , S u (ω) will have an extremum at ω 2 = z + , provided Otherwise, it will be monotonically decreasing from the value of S u (0) = A/C 4 to S u (∞) = 0. When the condition (24) holds, S u (ω) will have its single extremum at ω 2 0 = z + , and, due to continuity, this extremum will be a maximum.
In order to better understand the structure of stochastic oscillations around the dominant spectral frequency ω 2 0 = A 2 + B(BC 4 + 2AC 2 − AD 2 ) − A /B, we will use the notion of coherence that can be defined as follows. We consider an interval of frequencies [ω 1 , ω 2 ] around the dominant frequency ω 0 , i.e., ω 1 < ω 0 < ω 2 , and compute a quantity (cf. expressions in Eqs. (2.8) and (2.9) in Ref. 89) and a factor of 2 in front of the integral is introduced to account for the fact that we will only be considering non-negative frequencies 0 ≤ ω < ∞, rather than the full range −∞ < ω < ∞. Evaluating the same integral over the entire real line yields and the coherence c can then defined as 89 which is a non-negative quantity satisfying c = 0 for ω 1 = ω 2 = ω 0 and c = 1 for ω 1 = −∞ and ω 2 = ∞. 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 * , 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
While stochastic oscillations in the neighborhood of the stable steady state E * provide some understanding of how small noise may affect otherwise decaying oscillations, an interesting and important question is whether, under small amount of noise, it would be possible for some stochastic trajectories to leave the neighborhood of E * and reach the separatrix, from which they could escape further and reach the right branch of u-nullcline, associated with an orgasm. We will, therefore, interpret any trajectory, for which the values of physiological arousal during time evolution exceed the value of u = u B , as an orgasmic trajectory. To explore the dynamics of such stochastic escape, we rescale noise parameters and change variables to In these new variables, the Langevin model (13) is transformed into a system where is the deterministic part of the model in new variables and W t , is a two-dimensional Wiener process. We are interested in solutions of this equation that will leave some region R ∈ R 2 enclosing the steady state E * after a certain time T. Since this steady state is deterministically stable, in the case when µ = 0, solutions would never leave its neighborhood. However, for small noise µ > 0, some of them will indeed escape, but the probability of this escape becomes smaller and smaller as µ → 0. Large deviation theory (LDT) 90 provides a quantitative characterization of this decay of probabilities. If we introduce an action functional, also known as the rate function, as with Lagrangian L(x,ẋ) being defined as follows, then for any chosen vector function x 0 (t), for small enough δ > 0, the probability of observing a sample path of (26) close to this x 0 (t) can be estimated as Here, denotes logarithmic equivalence: We can now define quasi-potential V(x 1 , x 2 ) 91-93 that characterizes long-term behavior of the system as where C x 1 ,x 2 = x ∈ AC([0, T], R 2 ), x(0) = x 1 , x(T) = x 2 is the space of absolutely continuous vector functions. Since our model has a single stable steady state E * , one can write probability density ρ(x) associated with the invariant density of (26) in the form of WKB approximation, 90 where This allows us to estimate the expected time of escape τ esc from the basin of attraction of E * , denoted B(E * ), as B(E * ), and it produces the smallest action S T [x * (t)]; i.e.,

Chaos
We note that the minimization problem defining the most probable escape path (31) directly corresponds to Hamilton's principle δS T [x]/δx = 0. Using this analogy, one can define "conjugate momenta" Minimization in Eq. (31) is then equivalent to finding solutions to Hamilton's equations of motion, or "instanton equations," 78,91 which have the explicit forṁ with the corresponding equation for the actioṅ Conjugate momenta can be interpreted as some quantitative measure of the deviation of most probable escape paths from the deterministic dynamics. Steady state E * of the original model corresponds to the steady state ( E * , 0) of the extended system (33), but its stability changes: while originally E * was a stable focus, steady state ( E * , 0) of the extended system (33) is a saddle. With separatrix playing the role of quasi-threshold, 59 MPEP can be viewed as solutions of the extended system that originate in the small neighborhood of the saddle ( E * , 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, 95 in which a large number of initial conditions are chosen to be uniformly distributed on a small circle around the steady state ( E * , 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), [91][92][93] 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 91 and the adaptive minimum action method, 96,97 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, 98 a Ludwig spruce budworm model, 78 a predator-prey model with migration, 91 and a Maier-Stein model 99 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 * , 0) to points on a separatrix.
While path-based methods, such as gMAM and an action plot, allow one to find MPEP, they have certain limitations, which include a bias introduced by the initial path, potential inexactness of a minimizer due to slow convergence, and the possibility that a found local action minimizer might not be a global minimizer. 100 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) 101 and an ordered line integral method (OLIM) 100,102 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, 103 the Zeldovich-Semenov model for an exothermic reaction in a continuously stirred tank reactor, 104 and the FitzHugh-Nagumo model of neural excitation, 59 among others. Once the quasi-potential U E * is found, MPEP that goes from the attractor E * to a given point x can be found by integrating the equation 100 from the starting point y(0) = x back to E * , and normalization is introduced to avoid slowdown of convergence near the steady state E * .
To compute quasi-potential and MPEP for the system (26), we considered a region [−1.5, 7] × [0, 11] that was divided by a rectangular mesh 2000 × 2000, over which quasi-potential was computed in MATLAB with a modification of the OUM code developed earlier for the Higgins model. 105 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 nonmonotonic, 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, 90 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 . . 20, 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. 95,106 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; 95 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 φ = 0.0014646947 to set up initial conditions and then computed MPEP by following the Hamiltonian system (33), as shown in Fig. 9. Having E * 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 quasipotential 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. [107][108][109] While the reasons for this can be both psychogenic and organic, it Chaos ARTICLE scitation.org/journal/cha has been shown that with damage to the spinal cord, there may be greater reliance on psychogenic stimulation for achieving arousal and orgasm, 110,111 and while the distinction between reflexive and psychogenic erection in able-bodied males may be less clear, 112 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. 113,114 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. 74 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. 115 A more recent survey 113 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, 116 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," 117 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. [118][119][120] 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, 121 chaotic synchronization, 122 and chimera states. 123 In this context, it has been suggested that periodic physiological Chaos ARTICLE scitation.org/journal/cha stimulation, which we modeled as a constant term due to its highfrequency nature when considered on the timescale of a sexual response cycle, can actually result in the entrainment of oscillations among groups of neurons, 118 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. 124,125 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.