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 sexuality4 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–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 arousal14–16 and an observation that arousal was continuing to increase during the plateau phase.17 Levin18,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
(a) Surface of steady states of the system (1). (b) Projection onto the – parameter plane illustrating a cusp, as well as fold curves separating parameter regions with one and two stable steady states. (c) Sample curves illustrating monotonic behavior (black, ) and hysteresis (red, ), with curves of stable steady states shown in solid and unstable in dashed.
(a) Surface of steady states of the system (1). (b) Projection onto the – parameter plane illustrating a cusp, as well as fold curves separating parameter regions with one and two stable steady states. (c) Sample curves illustrating monotonic behavior (black, ) and hysteresis (red, ), with curves of stable steady states shown in solid and unstable in dashed.
Hubey29,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 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.
Numerical solution of Duffing equation (2) with , , and . (a) Periodic solution for . (b) Period-two solution for . (c) Period-four solution for . (d) Chaotic solution for .
Numerical solution of Duffing equation (2) with , , and . (a) Periodic solution for . (b) Period-two solution for . (c) Period-four solution for . (d) Chaotic solution for .
Levin18,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
(a) Nullclines of the system (6)–(8). (b) -nullcline with (blue) or (red).
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 with some very high frequency , we rather represent it as a constant over the entire duration of the sex cycle.
The term 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 amygdala38 and temporal lobe,39 as well as with deactivation of frontal cortical regions22,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 as stated in that equation or, alternatively, as a cubic function of the form , 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 . 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).
IV. STABILITY ANALYSIS AND NUMERICAL SIMULATIONS
When is lying on the middle branch, i.e., , the condition of uniqueness of the steady state guarantees that . Then, the steady state is stable for , unstable for , and undergoes a supercritical Hopf bifurcation when . We note that since for all values of , this means that is always positive, and hence, 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 , 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 will regain its stability. Finally, on the right branch , the steady state is stable. The points and are always stable, and they are either stable nodes if or stable foci otherwise. In light of condition (10) that ensures uniqueness of the steady state and the fact that , in most cases, the points and will be stable foci.
Bifurcation diagram for the steady state of the model (6)–(8) with , . (a) , . (b) , . Colors denote a stable node (black), a stable focus (blue), an unstable focus (red), and an unstable node (yellow).
Bifurcation diagram for the steady state of the model (6)–(8) with , . (a) , . (b) , . Colors denote a stable node (black), a stable focus (blue), an unstable focus (red), and an unstable node (yellow).
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 , as well as stability type of the steady state . Figure 5(a) illustrates the sequence of transitions shown in Fig. 4 when the level of physiological stimulation is fixed, while the level of psychological stimulation is increasing. This corresponds to moving -nullcline upward, while -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 does not alter the shape of the bifurcation curve but results in higher values of and stability changes occurring at slightly higher values of . Figure 5(b) illustrates a similar sequence of transitions, though starting all the way on the right branch of the -nullcline and moving to the left along this nullcline for increasing values of .
One-dimensional bifurcation diagrams showing a steady-state value of depending on parameters. Colors correspond to those in Fig. 4 and denote an unstable node (yellow), an unstable focus (red), a stable focus (blue), and a stable node (black). Parameter values are , , . (a) . (b) , .
One-dimensional bifurcation diagrams showing a steady-state value of depending on parameters. Colors correspond to those in Fig. 4 and denote an unstable node (yellow), an unstable focus (red), a stable focus (blue), and a stable node (black). Parameter values are , , . (a) . (b) , .
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 -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 without ever reaching the right branch of the -nullcline. We define as the critical point in the sense that if the level of physiological arousal has exceeded this value before returning to , this will be interpreted as an orgasm having happened. Numerically, separatrix can be found by either taking initial conditions on the right branch of -nullcline and integrating the system backward in time58 or by taking the single point 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
(a) – phase plane of the model (6)–(8). Dashed black curves indicate nullclines, blue and red curves are anorgasmic and orgasmic trajectories, respectively, and the green curve indicates the separatrix. (b) Numerical solution of the model (6)–(8), red and blue curves correspond to phase trajectories in plot (a), and the dashed black line denotes an orgasmic level of physiological arousal. (c) Close-up of plot (b), indicating a period of initial fast-growing excitation (E), followed by a plateau (P), where arousal continues to increase, but at a slower pace, leading to either an orgasm (O) followed by resolution (R) or resolution without an orgasm. Parameter values are , , , , , .
(a) – phase plane of the model (6)–(8). Dashed black curves indicate nullclines, blue and red curves are anorgasmic and orgasmic trajectories, respectively, and the green curve indicates the separatrix. (b) Numerical solution of the model (6)–(8), red and blue curves correspond to phase trajectories in plot (a), and the dashed black line denotes an orgasmic level of physiological arousal. (c) Close-up of plot (b), indicating a period of initial fast-growing excitation (E), followed by a plateau (P), where arousal continues to increase, but at a slower pace, leading to either an orgasm (O) followed by resolution (R) or resolution without an orgasm. Parameter values are , , , , , .
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 , 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
A. Stochastic oscillations near equilibrium
(a) Phase plane of the system showing confidence ellipses (20) for with (blue), (red), and (magenta). The green curve illustrates one stochastic realization of the system (13). Dashed black curves show nullclines of the deterministic model (6)–(8). (b) Comparison of numerical solutions of the deterministic model (black) and its stochastic counterpart (green), given by Eqs. (6)–(8) and (13), respectively, with initial conditions in the neighborhood of the steady state . Parameter values are , , , , , , and .
(a) Phase plane of the system showing confidence ellipses (20) for with (blue), (red), and (magenta). The green curve illustrates one stochastic realization of the system (13). Dashed black curves show nullclines of the deterministic model (6)–(8). (b) Comparison of numerical solutions of the deterministic model (black) and its stochastic counterpart (green), given by Eqs. (6)–(8) and (13), respectively, with initial conditions in the neighborhood of the steady state . Parameter values are , , , , , , and .
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 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 of relaxation of psychological arousal to its state of homeostasis makes the frequency distribution sharper around the dominant frequency. Reducing psychological stimulation or reducing the rate , at which psychological arousal increases with physiological stimulation, leads to a higher variance of stochastic oscillations around the deterministically stable steady state , 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.
(a) Power spectrum of from (23) with , , , , and (black), (blue), (red), (green). (b) Power spectrum of from (23) with , , , , and (black), (blue), (red), (green). Variance (c) and (d) and coherence (e) and (f) of stochastic oscillations around the steady state computed from (17) and (25), respectively. Parameter values are , , for (c) and (e) and , , for plots (d) and (f).
(a) Power spectrum of from (23) with , , , , and (black), (blue), (red), (green). (b) Power spectrum of from (23) with , , , , and (black), (blue), (red), (green). Variance (c) and (d) and coherence (e) and (f) of stochastic oscillations around the steady state computed from (17) and (25), respectively. Parameter values are , , for (c) and (e) and , , for plots (d) and (f).
B. Large deviations: An orgasm as an escape
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 , 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 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 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 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 gMAM92 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 model100 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 to points on a separatrix.
To compute quasi-potential and MPEP for the system (26), we considered a region that was divided by a rectangular mesh , 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 – 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 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.
(a) Logarithm of the quasi-potential computed using OUM. (b) Contour lines of quasi-potential. (c) Action plot. (d) Optimal path (MPEP) computed using quasi-potential (black), action plot (magenta), and gMAM (blue). In all plots, dashed black curves indicate nullclines, and the solid green curve denotes separatrix. The red square, denoted as in plots (b) and (d), indicates the location of optimal escape, characterized by the lowest value of quasi-potential along the separatrix. Parameter values are , , , , and .
(a) Logarithm of the quasi-potential computed using OUM. (b) Contour lines of quasi-potential. (c) Action plot. (d) Optimal path (MPEP) computed using quasi-potential (black), action plot (magenta), and gMAM (blue). In all plots, dashed black curves indicate nullclines, and the solid green curve denotes separatrix. The red square, denoted as in plots (b) and (d), indicates the location of optimal escape, characterized by the lowest value of quasi-potential along the separatrix. Parameter values are , , , , and .
For an action plot, we took 20 000 initial conditions chosen as , where , , , uniformly distributed around the circle of radius centered on the steady state, with initial conditions for conjugate momenta and action 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 , associated with the minimum of quasi-potential along the separatrix, differs from the overall minimum value of action by less than , and hence, this point 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 to set up initial conditions and then computed MPEP by following the Hamiltonian system (33), as shown in Fig. 9. Having and 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 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 frequency5) and different societal groups.116 A more recent survey114 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 Johnson7 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 Basson8 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.