Sex, ducks, and rock “n” roll: Mathematical model of sexual response

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

Significant progress in this direction was made in the twentieth century, starting with the works of Freud on origins of sexuality⁴ 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,⁷ 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,”¹² and the dual control model.¹³ Some of the criticisms of the original Masters–Johnson raised in these later works included lack of contribution from the psychological component of the arousal^14–16 and an observation that arousal was continuing to increase during the plateau phase.¹⁷ Levin^18,19 provides a nice overview and comparison of these and other models of a human sexual response.

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

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

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 two-dimensional systems. In the former case, once the forcing frequency is fixed, this will result in some fixed amplitude that will not be changing with time. To move along the amplitude response curve, one would actually have to change a forcing frequency, and then for each value of the frequency, there would be an associated value of the amplitude of oscillations. In contrast, closed phase trajectories in the phase plane do represent periodic solutions once system parameters permit existence of such solutions for certain parameter values. There is no cusp catastrophe associated with the tip of the amplitude response curve, because even though at that point, a transition between branches occur, in each case, the system would exhibit oscillations with only one value of the amplitude, and there is no change in the number of steady states of the model. Moreover, whereas the Duffing oscillator can exhibit periodic, quasiperiodic, and chaotic behavior, as shown in Fig. 2, none of its solutions has the form that would qualitatively resemble temporal dynamics of the Masters–Johnson sexual response cycle. Furthermore, the Masters–Johnson sexual response cycle is concerned with a single trajectory showing growth of arousal culminating in an orgasm, while any point on the amplitude response curve rather corresponds to sustained periodic oscillations of a certain fixed amplitude. All this suggests that while there may be some scope in modeling the dynamics of sexual arousal as a forced two-dimensional system, associating hysteresis on the amplitude response curve with a cusp catastrophe is not justified.

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

With the average duration of a sexual intercourse being estimated to be around 5.4 min,³⁶ which includes between 100 and 500 thrusts,³⁷ instead of considering physical stimulation as a periodic forcing $E ucos⁡(ω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³⁸ and temporal lobe,³⁹ as well as with deactivation of frontal cortical regions^22,40,41 and orbitofrontal cortex.²² This suggests that reaching an orgasm can be characterized as psychologically “letting go,”⁴² 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,⁴⁵ 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 $(−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).

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( J 0)=ϵ [ a − b f ′ ( u ∗ ) ]>0$⁠. Then, the steady state $E ∗$ is stable for $f ′( u ∗)−bϵ<0$⁠, unstable for $f ′( u ∗)−bϵ>0$⁠, and undergoes a supercritical Hopf bifurcation when $f ′( u ∗)−bϵ=0$⁠. We note that since $a−b f ′(u)>0$ for all values of $u$⁠, this means that $det( J 0)=ϵ [ a − b f ′ ( 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 Fig. 5, we present one-dimensional bifurcation diagrams obtained by fixing all of the parameters except one, and then for each value of bifurcation parameter, we note the corresponding steady-state value of psychological arousal $v ∗$⁠, as well as stability type of 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$⁠. Figure 5(b) illustrates a similar sequence of transitions, though starting all the way on the right branch of the $u$-nullcline and moving to the left along this nullcline for increasing values of $a$⁠.

Figure 6(a) illustrates the phase plane of the system (6)–(8). Similarly to the well-known case of a FitzHugh–Nagumo model, the left and right branches of the $u$-nullcline are attracting, while the middle branch is repelling, and there exists a separatrix shown in green such that initial conditions that lie in the phase space to the left of this separatrix return to the steady state $E ∗$ 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⁵⁸ or by taking the single point $(u,v)=( u B,f( u B)+ E u)$ as an initial condition and integrating the system backward in time. Importantly, in the same way as with the FitzHugh–Nagumo model, there is no clearly defined threshold in the model, and hence, separatrix plays the role of a quasi-threshold manifold.⁵⁹ 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

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.⁷⁵ 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 

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.

One method of computing MPEP is the method of an action plot,⁹⁶ 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),^92–94 in which the fact that the action $S( x)$ is independent of path parameterization allows one to replace that action by an action that is parameterized by normalized arc length instead of time. This allows one to overcome the numerical difficulty associated with the fact that since the escape trajectory starts at the steady state, time $T$ in the definition of action is infinite, whereas it would be replaced by a finite arc length once action is replaced by a geometric action. Several versions of the gMAM method have been proposed, including the simplified gMAM⁹² 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,⁵⁹ a system with overdamped double-well potential,⁹⁹ a Ludwig spruce budworm model,⁷⁹ a predator–prey model with migration,⁹² and a Maier–Stein model¹⁰⁰ 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.

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.¹⁰⁶ The resulting quasi-potential and its contour lines over the $u ~$– $v ~$ plane are shown in Figs. 9(a) and 9(b). The shape of quasi-potential can be described as a deep well around the steady state, with a plateau-like area further away. Since we are interested in the MPEP, where stochastic trajectories escape a stable steady state by crossing a quasi-threshold associated with the separatrix, we investigated how the value of quasi-potential changes along this separatrix. This analysis reveals that the change is non-monotonic, and there is a point that we denote as $S 0$ that represents the minimum of quasi-potential along the separatrix. Since, according to (30), such a point would provide the shortest escape time and would also correspond to the optimal escape path,⁹¹ 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.