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

Blyuss, K. B.; Kyrychko, Y. N.

doi:10.1063/5.0143190

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

I. INTRODUCTION

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

Significant progress in this direction was made in the twentieth century, starting with the works of Freud on origins of sexuality⁴ 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 $\to$ plateau $\to$ orgasm $\to$ 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.

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

\dot{x} = \frac{d V}{d x},

(1)

with

V (x) = x^{4} - a x - \frac{b}{2} x^{2} .

(2)

Steady states of this model satisfy a cubic equation

x^{3} = a + b x,

(3)

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.

FIG. 1.

View large Download slide

(a) Surface of steady states of the system (1). (b) Projection onto the $a$ – $b$ 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, $b = - 1$ ⁠) and hysteresis (red, $b = 1$ ⁠), with curves of stable steady states shown in solid and unstable in dashed.

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,

\ddot{θ} + α \dot{θ} + β θ + γ θ^{3} = f \cos (Ω t),

(4)

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³⁴

| θ |^{2} [{(α Ω)}^{2} + {(Ω^{2} - β - \frac{3}{4} γ | θ |^{2})}^{2}] = f^{2} .

(5)

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

FIG. 2.

View large Download slide

Numerical solution of Duffing equation (2) with $α = 0.3$ ⁠, $β = - 1$ ⁠, and $γ = 1$ ⁠. (a) Periodic solution for $f = 0.25$ ⁠. (b) Period-two solution for $f = 0.32$ ⁠. (c) Period-four solution for $f = 0.35$ ⁠. (d) Chaotic solution for $f = 0.38$ ⁠.

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³⁵ that a sexual response is a complex process with multiple contributions, we model sexual arousal by two variables,

u (t)

and

v (t)

⁠, which are taken to represent current levels of physiological and psychological/cognitive arousal, respectively. In doing so, we address some of the above-mentioned criticisms of the Masters–Johnson model in terms of not appropriately accounting for the psychological component of the arousal.^14–16 The equation for the physiological component of arousal

u

is taken to be of the form

\dot{u} = f (u) - v + E_{u},

(6)

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 (\cdot)

is chosen to be of the following form that is inspired by earlier works on excitable systems:

f (u) = {\begin{cases} f_{L} (u), u \leq u_{A}, \\ f_{M} (u), u_{A} \leq u \leq u_{B}, \\ f_{R} (u), u \geq u_{B}, \end{cases}

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

f_{L} (u_{A}) = f_{M} (u_{A})

⁠,

f_{M} (u_{B}) = f_{R} (u_{B})

⁠, and

f_{L}^{'} (u_{A}) = f_{M}^{'} (u_{A}) = f_{M}^{'} (u_{B}) = f_{R}^{'} (u_{B}) = 0

⁠, they will still provide a continuously differentiable right-hand side of Eq. . We will choose these functions to be

f (u) = {\begin{cases} f_{L} (u) = f_{M} (u) = f_{1} (u) = a_{0} + a_{1} u + a_{2} u^{2} + a_{3} u^{3}, \\ f_{R} (u) = f_{2} (u) = f_{M} (u_{B}) - {(u - u_{B})}^{5} / 5 \\ = b_{0} + b_{1} u + b_{2} u^{2} + b_{3} u^{3} + b_{4} u^{4} + b_{5} u^{5} . \end{cases}

(7)

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.

FIG. 3.

View large Download slide

(a) Nullclines of the system (6)–(8). (b) $u$ -nullcline with $f (u) + E_{u} = v$ (blue) or $f (u) + E_{u} + v - v^{3} / 3 = 0$ (red).

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_{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³⁸ 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).

Dynamics of a psychological component of arousal is modeled by the following equation:

\dot{v} = ϵ [(E_{v} - E_{v_{0}}) + a u - b v],

(8)

where

E_{v 0}

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

a u

⁠. Justification for the last term

(- b v)

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–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. ], and the B-process that is responsible for returning to the state of psychological homeostasis. This theory has already been effectively used to study a number of psychological processes, such as drug addiction and withdrawal,^52–55 as well as sleep.^56,57 Finally, whereas physiological responses are rather quick and are controlled by the autonomic nervous system, processing of psychological stimuli is (much) slower; hence, we introduce a scaling parameter

0 < ϵ ≪ 1

⁠. In both Eqs. and , parameters

E_{v 0}

⁠,

a

⁠, and

b

are assumed to be positive, while external stimulations

E_{u}

and

E_{v}

are non-negative.

IV. STABILITY ANALYSIS AND NUMERICAL SIMULATIONS

In order to understand the structure of steady states of the model , it is instructive to explore its nullclines

\begin{array}{l} \dot{u} = 0 ⟹ v = f (u) + E_{u}, \\ \dot{v} = 0 ⟹ v = \frac{1}{b} [(E_{v} - E_{v_{0}}) + a u] . \end{array}

(9)

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:

\frac{a}{b} > f^{'} (u), \forall u \in [u_{A}, u_{B}],

(10)

which guarantees that the model always has a single steady state

E^{*} = (u^{*}, v^{*})

⁠. Introducing auxiliary parameters

k = a / b

and

d = (E_{v 0} - E_{v}) / a

⁠, we then have the following classification:

\begin{aligned} k (u_{A} - d) > f (u_{A}) + E_{u}, E^{*} lies on the left branch, \\ k (u_{A} - d) = f (u_{A}) + E_{u}, E^{*} = (u_{A}, f (u_{A}) + E_{u}), \\ k (u_{A} - d) < f (u_{A}) + E_{u}, k (u_{B} - d) > f (u_{B}) + E_{u}, \\ E^{*} lies on the middle branch, \\ k (u_{B} - d) = f (u_{B}) + E_{u}, E^{*} = (u_{B}, f (u_{B}) + E_{u}), \\ k (u_{A} - d) < f (u_{A}) + E_{u}, k (u_{B} - d) < f (u_{B}) + E_{u}, \\ E^{*} lies on the right branch . \end{aligned}

For the particular choice of the cubic function

f

given in , the locations of its two critical points are given by

u_{A, B} = \frac{- a_{2} \mp \sqrt{a_{2}^{2} - 3 a_{1} a_{3}}}{3 a_{3}} .

Furthermore, since this function is cubic, it will have an inflection point

u_{i n f l} = - \frac{a_{2}}{3 a_{3}},

where

f^{″} (u_{i n f l}) = 0

and

f^{'} (u_{i n f l})

reaches its maximum; therefore, the condition for uniqueness of the steady state turns into

\frac{a}{b} > f^{'} (u_{i n f l}) = a_{1} - \frac{a_{2}^{2}}{3 a_{3}} .

(11)

Jacobian of linearization near the steady state

E^{*}

is given by

J_{0} = (\begin{array}{cc} f^{'} (u^{*}) & - 1 \\ a ϵ & - b ϵ \end{array}),

(12)

with

Tr (J_{0}) = f^{'} (u^{*}) - b ϵ

and

det (J_{0}) = ϵ [a - b f^{'} (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

E^{*}

is stable. If

a

is sufficiently small so that

{[f^{'} (u^{*}) + b ϵ]}^{2} - 4 ϵ a \geq 0

⁠,

E^{*}

is a stable node, otherwise, it is a stable focus, and at

a = {[f^{'} (u^{*}) + b ϵ]}^{2} / 4 ϵ

⁠, there are two coinciding characteristic roots with values

λ_{1} = λ_{2} = [f^{'} (u^{*}) - b ϵ] / 2 < 0

⁠.

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} ϵ \geq 4 a$ 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

f (u) = {\begin{cases} f_{1} (u) & = \frac{1}{2} (u + \frac{3}{2}) - \frac{3}{2} - \frac{1}{3} {[\frac{1}{2} (u + \frac{3}{2}) - \frac{5}{2}]}^{3} \\ = \frac{167}{192} - \frac{33}{32} u + \frac{7}{16} u^{2} - \frac{1}{24} u^{3}, u \leq u_{B}, \\ f_{2} (u) & = \frac{3}{2} - {(u - \frac{11}{2})}^{5} = \frac{161099}{32} - \frac{73205}{16} u \\ + \frac{6655}{4} u^{2} - \frac{605}{2} u^{3} + \frac{55}{2} u^{4} - u^{5}, u \geq u_{B}, \end{cases}

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 - b f^{'} (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.

FIG. 4.

View large Download slide

Bifurcation diagram for the steady state $E^{*}$ of the model (6)–(8) with $ϵ = 0.05$ ⁠, $E_{v_{0}} = 7 / 4$ ⁠. (a) $a = 0.5$ ⁠, $b = 0.1$ ⁠. (b) $E_{u} = 1 / 6$ ⁠, $E_{v} = 1.05$ ⁠. 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 $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$ ⁠.

FIG. 5.

View large Download slide

One-dimensional bifurcation diagrams showing a steady-state value of $v^{*}$ 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 $ϵ = 0.05$ ⁠, $b = 0.1$ ⁠, $E_{v_{0}} = 7 / 4$ ⁠. (a) $a = 0.5$ ⁠. (b) $E_{u} = 1 / 6$ ⁠, $E_{v} = 1.05$ ⁠.

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

FIG. 6.

View large Download slide

(a) $u$ – $v$ 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 $ϵ = 0.05$ ⁠, $a = 0.5$ ⁠, $b = 0.1$ ⁠, $E_{u} = 1 / 6$ ⁠, $E_{v} = 1.05$ ⁠, $E_{v_{0}} = 7 / 4$ ⁠.

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

In light of these observations, we modify the model by writing it as a system of Langevin equations,

\begin{matrix} \dot{u} = f (u) - v + E_{u} + \sqrt{2 D_{1}} η_{1} (t), \\ \dot{v} = ϵ [(E_{v} - E_{v 0}) + a u - b v] + \sqrt{ϵ} \sqrt{2 D_{2}} η_{2} (t), \end{matrix}

(13)

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^{'}) ⟩ = δ_{i j} δ (t - t^{'})

for

i, j = 1, 2

⁠. An extra factor

\sqrt{ϵ}

in the noise term in the second equation of stochastic model has been introduced to ensure noises act on their respective timescales.^79,80

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 near this deterministically stable steady state and obtain a system of equations for stochastic fluctuations

ξ_{1, 2} (t)

⁠,

\begin{matrix} \frac{d}{d t} ξ = J_{0} ξ + B η, ξ = (\begin{matrix} ξ_{1} \\ ξ_{2} \end{matrix}), η = (\begin{matrix} η_{1} \\ η_{2} \end{matrix}), \\ B = (\begin{array}{cc} \sqrt{2 D_{1}} & 0 \\ 0 & \sqrt{2 D_{2} ϵ} \end{array}), \end{matrix}

(14)

where

J_{0}

is the Jacobian 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 describes a two-dimensional Ornstein–Uhlenbeck process.^81–83 In this case, the covariance matrix

Σ

for fluctuations is constant and is given by the solution of the Lyapunov equation,⁸³ known in the context of control problems as the continuous-time algebraic Riccati equation,⁸⁴

J_{0} Σ + Σ J_{0}^{T} + 2 D = 0,

(15)

where

D

is the diffusion matrix,

D = \frac{1}{2} B B^{T} = (\begin{array}{cc} D_{1} & 0 \\ 0 & D_{2} ϵ \end{array}) .

The solution of this equation can be found in a closed form,^81,83

\begin{aligned} Σ & = ⟨ ξ (t) ξ^{T} (t^{'}) ⟩ = - \frac{det (J_{0}) D + [J_{0} - Tr (J_{0}) I] D {[J_{0} - Tr (J_{0}) I]}^{T}}{Tr (J_{0}) det (J_{0})} \\ = - \frac{1}{f^{'} (u^{*}) - b ϵ} [(\begin{array}{cc} D_{1} & 0 \\ 0 & D_{2} ϵ \end{array}) + \frac{1}{a - b f^{'} (u^{*})} \\ \times (\begin{array}{cc} D_{1} b^{2} ϵ + D_{2} & D_{1} a b ϵ + D_{2} f^{'} (u^{*}) \\ D_{1} a b ϵ + D_{2} f^{'} (u^{*}) & D_{1} a^{2} ϵ + D_{2} {[f^{'} (u^{*})]}^{2} \end{array})] . \end{aligned}

(16)

This gives the variance of stochastic fluctuations in the level of physiological arousal as

σ_{u}^{2} = Σ_{u u} = - \frac{1}{f^{'} (u^{*}) - b ϵ} [D_{1} + \frac{D_{1} b^{2} ϵ + D_{2}}{a - b f^{'} (u^{*})}] .

(17)

The Fokker–Planck equation for the probability density

ρ (ξ, t)

associated with system has the form

\frac{\partial ρ}{\partial t} = - \nabla (ρ J_{0}) + \sum_{i, j} \frac{\partial^{2}}{\partial ξ_{i} \partial ξ_{j}} (D_{i j} ρ) .

(18)

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,^81,83

ρ (ξ, t) = \frac{1}{2 π \sqrt{det (J_{0})}} \exp (- \frac{1}{2} ξ^{T} Σ^{- 1} ξ) .

(19)

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).^85–87 Using SSF, we can find the so-called confidence ellipse

{(x - x_{0})}^{T} Σ^{- 1} (x - x_{0}) = 2 k^{2}, x = (\begin{array}{l} u \\ v \end{array}), x_{0} = (\begin{array}{l} u^{*} \\ v^{*} \end{array}),

(20)

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⁸⁸^,

\frac{ζ_{1}^{2}}{μ_{1}} + \frac{ζ_{2}^{2}}{μ_{2}} = k^{2},

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 , illustrating the phenomenon of stochastic amplification,^89,90 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.⁷⁵

FIG. 7.

View large Download slide

(a) Phase plane of the system showing confidence ellipses (20) for $p = 0.95$ with $D_{1} = D_{2} = 10^{- 4}$ (blue), $D_{1} = D_{2} = 10^{- 3}$ (red), and $D_{1} = D_{2} = 10^{- 2}$ (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 $E^{*}$ ⁠. Parameter values are $ϵ = 0.05$ ⁠, $a = 0.5$ ⁠, $b = 0.1$ ⁠, $E_{u} = 1 / 6$ ⁠, $E_{v} = 1.05$ ⁠, $E_{v_{0}} = 7 / 4$ ⁠, and $D_{1} = D_{2} = 10^{- 4}$ ⁠.

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 , which yields

\begin{matrix} - i ω \tilde{ξ} (ω) = J_{0} \tilde{ξ} (ω) + B \tilde{η} (ω) \\ ⟹ \tilde{ξ} (ω) = - {(J_{0} + i ω I)}^{- 1} B \tilde{η} (ω) . \end{matrix}

This then gives the power spectrum matrix as

\begin{aligned} S (ω) & = ⟨ \tilde{ξ} (ω) \tilde{ξ} {(ω)}^{†} ⟩ = {(J_{0} + i ω I)}^{- 1} B B^{T} {(J_{0}^{T} - i ω I)}^{- 1} \\ = \frac{2}{{[ω^{2} - det (J_{0})]}^{2} + {[Tr (J_{0})]}^{2} ω^{2}} \\ \times (\begin{array}{cc} D_{1} (a_{22}^{2} + ω^{2}) + D_{2} ϵ a_{12}^{2} & - D_{1} a_{21} \\ - D_{2} ϵ a_{12} & D_{1} a_{21}^{2} + D_{2} ϵ (a_{11}^{2} + ω^{2}) \end{array}), \end{aligned}

(21)

where

†

denotes a Hermitian conjugate, and, for notational convenience, we have denoted as

a_{i j}

coefficients of the Jacobian matrix

J_{0} = (\begin{array}{cc} a_{11} & a_{12} \\ a_{21} & a_{22} \end{array}) .

(22)

For stochastic fluctuations in physiological arousal, we then have

\begin{aligned} S_{u} (ω) & = 2 \frac{D_{1} b^{2} ϵ^{2} + D_{2} ϵ + D_{1} z}{{(z - ϵ [a - b f^{'} (u^{*})])}^{2} + {[f^{'} (u^{*}) - b ϵ]}^{2} z} \\ = \frac{A + B z}{{(z - C^{2})}^{2} + D^{2} z}, z = ω^{2}, \end{aligned}

(23)

with

\begin{matrix} A = 2 ϵ (D_{1} b^{2} ϵ + D_{2}), C^{2} = det (J_{0}) = ϵ [a - b f^{'} (u^{*})], \\ B = 2 D_{1}, D = Tr (J_{0}) = f^{'} (u^{*}) - b ϵ, \end{matrix}

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

z_{\pm} = - \frac{A \pm \sqrt{A^{2} + B (B C^{4} + 2 A C^{2} - A D^{2})}}{B},

and since we require

z > 0

due to the relation

z = ω^{2}

⁠,

S_{u} (ω)

will have an extremum at

ω^{2} = z_{+}

⁠, provided

B C^{4} + 2 A C^{2} - A D^{2} > 0.

(24)

Otherwise, it will be monotonically decreasing from the value of

S_{u} (0) = A / C^{4}

to

S_{u} (\infty) = 0

⁠. When the condition holds,

S_{u} (ω)

will have its single extremum at

ω_{0}^{2} = 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

ω_{0}^{2} = [\sqrt{A^{2} + B (B C^{4} + 2 A C^{2} - A D^{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. 90)

P_{ω_{1}, ω_{2}} = 2 \int_{ω_{1}}^{ω_{2}} S_{u} (ω) d ω = F (ω_{2}) - F (ω_{1}),

where

\begin{aligned} F (ω) & = \frac{\sqrt{2}}{D \sqrt{D^{2} - 4 C^{2}}} (\frac{2 A + B D (\sqrt{D^{2} - 4 C^{2}} - D) + 2 B C^{2}}{\sqrt{D (D - \sqrt{D^{2} - 4 C^{2}}) - 2 C^{2}}} \\ \times \arctan \frac{ω \sqrt{2}}{\sqrt{D (D - \sqrt{D^{2} - 4 C^{2}}) - 2 C^{2}}} \\ - \frac{2 A - B D (\sqrt{D^{2} - 4 C^{2}} + D) + 2 B C^{2}}{\sqrt{D (D + \sqrt{D^{2} - 4 C^{2}}) - 2 C^{2}}} \\ \times \arctan \frac{ω \sqrt{2}}{\sqrt{D (D + \sqrt{D^{2} - 4 C^{2}}) - C^{2}}}), \end{aligned}

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 \leq ω < \infty

⁠, rather than the full range

- \infty < ω < \infty

⁠. Evaluating the same integral over the entire real line yields

P_{0} = \int_{- \infty}^{\infty} S_{u} (ω) d ω = \frac{π (A + B C^{2})}{D C^{2}},

and the coherence

c

can then defined as⁹⁰

c = \frac{P_{ω_{1}, ω_{2}}}{P_{0}},

(25)

which is a non-negative quantity satisfying

c = 0

for

ω_{1} = ω_{2} = ω_{0}

and

c = 1

for

ω_{1} = - \infty

and

ω_{2} = \infty

⁠.

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.

FIG. 8.

View large Download slide

(a) Power spectrum of $S_{u} (ω)$ from (23) with $ϵ = 0.05$ ⁠, $a = 0.5$ ⁠, $b = 0.1$ ⁠, $E_{u} = 0.3$ ⁠, and $E_{v} = 1$ (black), $E_{v} = 1.13$ (blue), $E_{v} = 1.16$ (red), $E_{v} = 1.2$ (green). (b) Power spectrum of $S_{u} (ω)$ from (23) with $ϵ = 0.05$ ⁠, $E_{u} = 1 / 6$ ⁠, $E_{v} = 1.05$ ⁠, $a = 2$ ⁠, and $b = 0.4$ (black), $b = 0.8$ (blue), $b = 1.2$ (red), $b = 1.6$ (green). Variance (c) and (d) and coherence (e) and (f) of stochastic oscillations around the steady state $E^{*}$ computed from (17) and (25), respectively. Parameter values are $ϵ = 0.05$ ⁠, $a = 0.5$ ⁠, $b = 0.1$ for (c) and (e) and $ϵ = 0.05$ ⁠, $E_{u} = 1 / 6$ ⁠, $E_{v} = 1.05$ for plots (d) and (f).

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

μ = 2 D_{1}, d = \sqrt{\frac{D_{2}}{D_{1}} ϵ}

and change variables to

x = (\begin{array}{l} \tilde{u} \\ \tilde{v} \end{array}), \tilde{u} = u, \tilde{v} = \frac{v}{d} .

In these new variables, the Langevin model is transformed into a system

d x = K (x) + \sqrt{μ} d W_{t},

(26)

where

K (x) = (\begin{matrix} f (\tilde{u}) - d \tilde{v} + E_{u} \\ \frac{ϵ}{d} [(E_{v} - E_{v 0}) + a \tilde{u} - b d \tilde{v}] \end{matrix})

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 \in 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

μ \to 0

⁠. Large deviation theory (LDT)⁹¹ provides a quantitative characterization of this decay of probabilities. If we introduce an action functional, also known as the rate function, as

S_{T} [x] = \int_{0}^{T} L (x, \dot{x}) d t,

(27)

with Lagrangian

L (x, \dot{x})

being defined as follows,

L (x, \dot{x}) = \frac{1}{2} ‖ \dot{x} - K (x) ‖^{2},

then for any chosen vector function

x_{0} (t)

⁠, for small enough

δ > 0

⁠, the probability of observing a sample path of close to this

x_{0} (t)

can be estimated as

P (sup_{0 \leq t \leq T} | x (t) - x_{0} (t) | < δ) ≍ \exp [- \frac{S_{T} [x_{0} (t)]}{μ}] .

Here,

≍

denotes logarithmic equivalence:

f (x)

is logarithmically equivalent to

\exp [α g (x)]

if

lim_{α \to \infty} \frac{1}{α} \log f (x) = g (x) .

We can now define quasi-potential

V (x_{1}, x_{2})

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

V (x_{1}, x_{2}) = inf_{T > 0} inf_{x \in C_{x_{1}, x_{2}}} S_{T} (x),

(28)

where

C_{x_{1}, x_{2}} = {x \in 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 in the form of WKB approximation,⁹¹

ρ (x) ≍ \exp [- \frac{U_{E^{*}} (x)}{μ}],

(29)

where

U_{E^{*}} (x) = V (E^{*}, x) .

This allows us to estimate the expected time of escape

τ_{e s c}

from the basin of attraction of

E^{*}

⁠, denoted

B (E^{*})

⁠, as

E [τ_{e s c}] ≍ \exp [min_{x \in \partial B (E^{*})} U_{E^{*}} (x) / μ] .

(30)

The most probable escape path (MPEP)

x^{*} (t)

⁠, also known as an instanton,^79,92,95 is effectively the least unlikely trajectory that escapes

B (E^{*})

⁠, and it produces the smallest action

S_{T} [x^{*} (t)]

⁠; i.e.,

x^{*} (t) = {argmin}_{x \in C} S_{T} [x],

(31)

where

C = {x \in AC ([0, T], R^{2}), x = x_{0}, x (T) \in B (E^{*})}

⁠.

We note that the minimization problem defining the most probable escape path directly corresponds to Hamilton’s principle

δ S_{T} [x] / δ x = 0

⁠. Using this analogy, one can define “conjugate momenta”

p = \frac{δ S_{T}}{δ x} = \frac{\partial L (x, \dot{x})}{\partial \dot{x}}

and perform the Fenchel–Legendre transform of the Lagrangian to obtain the Freidlin–Wentzell Hamiltonian,⁹¹^,

\begin{aligned} H (x, p) & = sup_{y} (⟨ p, y ⟩ - L (x, y)) \\ = \frac{1}{2} ⟨ p, p ⟩ + ⟨ K (x), p ⟩ = \frac{1}{2} p^{T} p + {[K (x)]}^{T} p . \end{aligned}

(32)

Minimization in Eq. is then equivalent to finding solutions to Hamilton’s equations of motion, or “instanton equations,”^79,92

\begin{matrix} \dot{x} = \nabla_{p} H (x, p) = K (x) + p, \\ \dot{p} = - \nabla_{x} H (x, p) = - {[\nabla K (x)]}^{T} p, \end{matrix}

which have the explicit form

\begin{matrix} \dot{\tilde{u}} = f (\tilde{u}) + E_{u} - d \tilde{v} + p_{u}, \\ \dot{\tilde{v}} = \frac{ϵ}{d} [(E_{v} - E_{v 0}) + a \tilde{u} - b d \tilde{v}] + p_{v}, \\ {\dot{p}}_{u} = - f^{'} (\tilde{u}) p_{u} - \frac{a ϵ}{d} p_{v}, \\ {\dot{p}}_{v} = d p_{u} + b ϵ p_{v}, \end{matrix}

(33)

with the corresponding equation for the action

\dot{S} = \frac{1}{2} (p_{u}^{2} + p_{v}^{2}) .

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

({\tilde{E}}^{*}, 0)

of the extended system , but its stability changes: while originally

E^{*}

was a stable focus, steady state

({\tilde{E}}^{*}, 0)

of the extended system is a saddle. With separatrix playing the role of quasi-threshold,⁵⁹ MPEP can be viewed as solutions of the extended system that originate in the small neighborhood of the saddle

({\tilde{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,⁹⁶ in which a large number of initial conditions are chosen to be uniformly distributed on a small circle around the steady state $({\tilde{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 $({\tilde{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.¹⁰¹ 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)¹⁰² and an ordered line integral method (OLIM)^101,103 have been proposed for numerically computing quasi-potential as a solution of the Hamilton–Jacobi equation associated with the Hamiltonian . OUM has been used to find quasi-potentials and to compute MPEP in the Higgins model for glycolysis,¹⁰⁴ the Zeldovich–Semenov model for an exothermic reaction in a continuously stirred tank reactor,¹⁰⁵ and the FitzHugh–Nagumo model of neural excitation,⁵⁹ 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^101–103

\dot{y} = - \frac{K (y) + \nabla U_{E^{*}} (y)}{‖ K (y) + \nabla U_{E^{*}} (y) ‖}, y (0) = x

(34)

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] \times [0, 11]$ that was divided by a rectangular mesh $2000 \times 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 $\tilde{u}$ – $\tilde{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.

FIG. 9.

View large Download slide

(a) Logarithm of the quasi-potential $U_{E^{*}} (x)$ 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 $S_{0}$ 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 $ϵ = 0.05$ ⁠, $a = 0.5$ ⁠, $b = 0.1$ ⁠, $E_{u} = 1 / 6$ ⁠, and $E_{v} = 1.05$ ⁠.

For an action plot, we took 20 000 initial conditions chosen as $x_{0 i} = E^{*} + δ x_{0 i}$ ⁠, where $δ x_{0 i} = {(R \cos ϕ_{i}, R \sin ϕ_{i})}^{T}$ ⁠, $R = 0.1$ ⁠, $i = 1 \dots 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.^96,107 Then, we solved the Hamiltonian system (33) for each of these initial conditions until the solution reached the separatrix, at which point the value of action was recorded. The resulting action plot is illustrated in Fig. 9(c), and it has the shape similar to that for other models;^96,104 namely, there is a minimum value of action, and there are several other initial conditions that yield values of action very close to the minimum. The value of action at the point $S_{0}$ ⁠, associated with the minimum of quasi-potential along the separatrix, differs from the overall minimum value of action by less than $0.04 %$ ⁠, and hence, this point $S_{0}$ can be chosen as the target point for computing MPEP in terms of minimizing action, as computed using an action plot. Hence, we chose the phase $ϕ = 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 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,¹¹³ 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,⁵ though in the majority of cases, both of these processes happen at the same time.⁷⁵ 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⁵) and different societal groups.¹¹⁶ A more recent survey¹¹⁴ 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,¹¹⁷ 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⁷ 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,”¹¹⁸ 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,¹²² chaotic synchronization,¹²³ and chimera states.¹²⁴ 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,¹¹⁹ thus suggesting the possibility of another feedback mechanism between physiological and psychological components of the arousal.

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

SUPPLEMENTARY MATERIAL

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

ACKNOWLEDGMENTS

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

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

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

DATA AVAILABILITY

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

REFERENCES

1.

J.

Keener

and

J.

Sneyd

,

Mathematical Physiology II: Systems Physiology

(

Springer

,

New York

,

2009

).

Google Scholar

Crossref

2.

J. T.

Ottesen

,

M. S.

Olufsen

, and

J. K.

Larsen

,

Applied Mathematical Models in Human Physiology

(

SIAM

,

Philadelphia, PA

,

2004

).

Google Scholar

Crossref

3.

S. J.

Chapman

,

A. C.

Fowler

,

R.

Hinch

, and

S. L.

Waters

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

4.

S.

Freud

,

Drei Abhandlungen Zur Sexualtheorie

(

F. Deuticke

,

Leipzig

,

1905

).

Google Scholar

5.

A. C.

Kinsey

,

W. R.

Pomeroy

, and

C. E.

Martin

,

Sexual Behavior in the Human Male

(

W.B. Saunders

,

Philadelphia, PA

,

1948

).

Google Scholar

6.

A. C.

Kinsey

,

W. R.

Pomeroy

,

C. E.

Martin

, and

P. H.

Gebhard

,

Sexual Behavior in the Human Female

(

W.B. Saunders

,

Philadelphia, PA

,

1953

).

Google Scholar

7.

W. H.

Masters

and

V. E.

Johnson

,

Human Sexual Response

(

Little, Brown and Co.

,

Boston, MA

,

1966

).

Google Scholar

8.

R.

Basson

,

J. Sex Marital Ther.

26

,

51

–

65

(

2000

).

https://doi.org/10.1080/009262300278641

Google Scholar

Crossref

PubMed

9.

R.

Basson

,

Br. J. Diabetes Vasc. Dis.

2

,

267

–

270

(

2002

).

https://doi.org/10.1177/14746514020020040501

Google Scholar

Crossref

10.

H. S.

Kaplan

,

J. Sex Marital Ther.

3

,

3

–

9

(

1977

).

https://doi.org/10.1080/00926237708405343

Google Scholar

Crossref

PubMed

11.

H. S.

Kaplan

,

Disorders of Sexual Desire

(

Simon & Schuster

,

New York

,

1979

).

Google Scholar

12.

R. J.

Levin

, “Human male sexuality: Appetite and arousal, desire and drive,” in Appetite-Neural and Behavioural Basis, edited by C. R. Legg and D. A. Booth (Oxford University Press, Oxford, 1994), pp. 127–164.

13.

J.

Bancroft

,

C. A.

Graham

,

E.

Janssen

, and

S. A.

Sanders

,

J. Sex Res.

46

,

121

–

142

(

2009

).

https://doi.org/10.1080/00224490902747222

Google Scholar

Crossref

PubMed

14.

J.

Bancroft

,

Human Sexuality and Its Problems

(

Churchill Livingstone Elsevier

,

London

,

2009

).

Google Scholar

15.

D.

Byrne

, “The imagery of sex,” in Handbook of Sexology, edited by J. Money and H. Musaph (Elsevier, Amsterdam, 1977), pp. 327–350.

16.

D.

Byrne

, “The study of sexual behavior as a multidisciplinary venture,” in Alternative Approaches to the Study of Sexual Behavior, edited by D. Byrne and H. Kelly (Erlbaum, Hillsdale, NJ, 1986), pp. 1–12.

17.

P. A.

Robinson

,

The Modernisation of Sex

(

Paul Elek Ltd.

,

London

,

1976

).

Google Scholar

18.

R. J.

Levin

, “The human sexual response cycle,” in The Textbook of Clinical Sexual Medicine, edited by W. W. IsHak (Springer, 2017), pp. 39-50.

19.

R. J.

Levin

,

Sex. Relat. Ther.

23

,

393

–

399

(

2008

).

https://doi.org/10.1080/14681990802488816

Google Scholar

Crossref

20.

N. J.

Wise

,

E.

Frangos

, and

B. R.

Komisaruk

,

J. Sex. Med.

14

,

1380

–

1391

(

2017

).

https://doi.org/10.1016/j.jsxm.2017.08.014

Google Scholar

Crossref

PubMed

21.

J. R.

Georgiadis

and

M. L.

Kringelbach

,

Prog. Neurobiol.

98

,

49

–

81

(

2012

).

https://doi.org/10.1016/j.pneurobio.2012.05.004

Google Scholar

Crossref

PubMed

22.

J. R.

Georgiadis

,

R.

Kortekaas

,

R.

Kuipers

,

A.

Nieuwenburg

,

J.

Pruim

,

A. A. T. S.

Reinders

, and

G.

Holstege

,

Eur. J. Neurosci.

24

,

3305

–

3316

(

2006

).

https://doi.org/10.1111/j.1460-9568.2006.05206.x

Google Scholar

Crossref

PubMed

23.

C.

Loehle

,

Ecol. Mod.

49

,

125

–

152

(

1989

).

https://doi.org/10.1016/0304-3800(89)90047-1

Google Scholar

Crossref

24.

D. D.

Jones

,

Simulation

29

,

1

–

15

(

1977

).

https://doi.org/10.1177/003754977702900102

Google Scholar

Crossref

25.

M. M.

Dodson

and

A.

Hallam

,

Am. Nat.

111

,

415

–

433

(

1977

).

https://doi.org/10.1086/283176

Google Scholar

Crossref

26.

D. A.

Rand

and

H. B.

Wilson

,

Proc. R. Soc. B

253

,

137

–

141

(

1993

).

https://doi.org/10.1098/rspb.1993.0093

Google Scholar

Crossref

27.

V.

Krane

,

Quest

44

,

72

–

87

(

1992

).

https://doi.org/10.1080/00336297.1992.10484042

Google Scholar

Crossref

28.

L.

Hardy

and

G.

Parfitt

,

Brit. J. Psychol.

82

,

163

–

178

(

1991

).

https://doi.org/10.1111/j.2044-8295.1991.tb02391.x

Google Scholar

Crossref

29.

H. M.

Hubey

, “Catastrophe theory and the human sexual response,” in 3rd International Symposium on Systems Research, Informatics and Cybernetics (Baden-Baden, Germany, 1991).

30.

H. M.

Hubey

, “Dynamics and geometry of the human sexual response,” in The Proceedings of the 2000 International Conference on Mathematics and Engineering Techniques in Medicine and Biological Sciences (METMBS2000) (Las Vegas, NV, 2000).

31.

R. J.

Levin

,

J. Sex. Med.

6

,

2376

–

2389

(

2009

).

https://doi.org/10.1111/j.1743-6109.2009.01350.x

Google Scholar

Crossref

PubMed

32.

C. A.

Isnard

and

E. C.

Zeeman

, “Some models from catastrophe theory in the social sciences,” in The Use of Models in the Social Sciences, edited by L. Collins (Routledge, New York, 1975).

33.

E. C.

Zeeman

,

Catastrophe Theory—Selected Papers 1972-1977

(

Addison-Wesley

,

Reading, MA

,

1977

).

Google Scholar

34.

D. W.

Jordan

and

P.

Smith

,

Nonlinear Ordinary Differential Equations—An Introduction for Scientists and Engineers

(

Oxford University Press

,

Oxford

,

2007

).

Google Scholar

35.

R. C.

Rosen

and

J. G.

Beck

,

Patterns of Sexual Arousal: Psychophysiological Processes and Clinical Applications

(

Guildford Press

,

New York

,

1988

).

Google Scholar

36.

E. W.

Corty

and

J. M.

Guardiani

,

J. Sex. Med.

5

,

1251

–

1256

(

2008

).

https://doi.org/10.1111/j.1743-6109.2008.00797.x

Google Scholar

Crossref

PubMed

37.

B. G.

Cooper

,

S. L.

Chin

,

R.

Xiao

,

K.

Buch

,

D.

Kim

, and

M. W.

Grinstaff

,

R. Soc. Open Sci.

5

,

180291

(

2018

).

https://doi.org/10.1098/rsos.180291

Google Scholar

Crossref

PubMed

38.

G.

Holstege

,

J. R.

Georgiadis

,

A. M. J.

Paans

,

L. C.

Meiners

,

F. H. C.

E. van der Graaf

, and

A. A. T.

Simone Reinders

,

J. Neurosci.

23

,

9185

–

9193

(

2003

).

https://doi.org/10.1523/JNEUROSCI.23-27-09185.2003

Google Scholar

Crossref

PubMed

39.

G.

Holstege

and

H. K.

Huynh

,

Horm. Behav.

59

,

702

–

770

(

2011

).

https://doi.org/10.1016/j.yhbeh.2011.02.008

Google Scholar

Crossref

PubMed

40.

J. R.

Georgiadis

,

A.

Reinders

,

A. M.

Paans

,

R.

Renken

, and

R.

Kortekaas

,

Hum. Brain Mapp.

30

,

3089

–

3101

(

2009

).

https://doi.org/10.1002/hbm.20733

Google Scholar

Crossref

PubMed

41.

J. R.

Georgiadis

,

M. J.

Farrell

,

R.

Boessen

et al.,

Neuroimage

50

,

208

–

216

(

2010

).

https://doi.org/10.1016/j.neuroimage.2009.12.034

Google Scholar

Crossref

PubMed

42.

K.

Sukel

, “Sex on the brain: Orgasms unlock altered consciousness,” New Scientist 2812 (2011); see http://www.newscientist.com/article/mg21028124.600-sex-on-the-brain-orgasms-unlock-alteredconsciousness.html?full=true.

43.

K. C.

Berridge

and

M. L.

Kringelback

,

Neuron

86

,

646

–

664

(

2015

).

https://doi.org/10.1016/j.neuron.2015.02.018

Google Scholar

Crossref

PubMed

44.

T.

Love

,

C.

Laier

,

M.

Brand

,

L.

Hatch

, and

R.

Hajela

,

Behav. Sci.

5

,

388

–

433

(

2015

).

https://doi.org/10.3390/bs5030388

Google Scholar

Crossref

45.

R. M.

Yerkes

and

J. D.

Dodson

,

J. Compar. Neurol. Psychol.

18

,

459

–

482

(

1908

).

https://doi.org/10.1002/cne.920180503

Google Scholar

Crossref

46.

G.

Patsis

,

H.

Sahli

,

W.

Verhels

, and

O.

De Troyer

, “Evaluation of attention levels in a tetris game using a brain computer interface,” in User Modeling, Adaptation, and Personalization, edited by S. Carberry, S. Weibelzahl, A. Micarelli, and G. Semeraro (Springer, Rome, 2013), pp. 127–138.

47.

S.

Arora

,

N.

Sevdalis

,

R.

Aggarwal

,

P.

Sirimanna

,

A.

Darzi

, and

R.

Kneebone

,

Urg. Endosc.

24

,

2588

–

2593

(

2010

).

https://doi.org/10.1007/s00464-010-1013-2

Google Scholar

Crossref

48.

C.

Stinson

and

D. A.

Bowman

,

IEEE Trans. Vis. Comput. Graph.

20

,

606

–

615

(

2014

).

https://doi.org/10.1109/TVCG.2014.23

Google Scholar

Crossref

PubMed

49.

R. L.

Solomon

and

J. D.

Corbit

,

Psychol. Rev.

81

,

119

–

145

(

1974

).

https://doi.org/10.1037/h0036128

Google Scholar

Crossref

PubMed

50.

R. L.

Solomon

and

J. D.

Corbit

,

Am. Econ. Rev.

68

,

12

–

24

(

1978

); available at https://www.jstor.org/stable/2951004.

Google Scholar

51.

R. L.

Solomon

,

Amer. Psychol.

35

,

691

–

712

(

1980

).

https://doi.org/10.1037/0003-066X.35.8.691

Google Scholar

Crossref

52.

T.

Chou

and

M. R.

D’Orsogna

,

Chaos

32

,

021102

(

2022

).

https://doi.org/10.1063/5.0082997

Google Scholar

Crossref

PubMed

53.

S. J.

Guastello

,

Beh. Sci.

29

,

258

–

262

(

1984

).

https://doi.org/10.1002/bs.3830290405

Google Scholar

Crossref

54.

G. F.

Koob

,

L.

Stinus

,

M.

Le Moal

, and

F. E.

Bloom

,

Neurosci. Biobehav. Rev.

13

,

135

–

140

(

1989

).

https://doi.org/10.1016/S0149-7634(89)80022-3

Google Scholar

Crossref

PubMed

55.

G. F.

Koob

,

S.

Barak Caine

,

L.

Parsons

,

A.

Markou

, and

F.

Weiss

,

Pharmacol. Biochem. Behav.

57

,

513

–

521

(

1997

).

https://doi.org/10.1016/S0091-3057(96)00438-8

Google Scholar

Crossref

PubMed

56.

D. M.

Edgar

,

W. C.

Dement

, and

C. A.

Fuller

,

J. Neurosci.

13

,

1065

–

1079

(

1993

).

https://doi.org/10.1523/JNEUROSCI.13-03-01065.1993

Google Scholar

Crossref

PubMed

57.

S. J.

Crowley

,

C.

Acebo

, and

M. A.

Carskadon

,

Sleep Med.

8

,

602

–

612

(

2007

).

https://doi.org/10.1016/j.sleep.2006.12.002

Google Scholar

Crossref

PubMed

58.

E. M.

Izhikevich

,

Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting

(

MIT Press

,

Cambridge, MA

,

2007

).

Google Scholar

59.

Z.

Chen

,

J.

Zhu

, and

X.

Liu

,

Proc. R. Soc. A

73

,

20170058

(

2017

).

https://doi.org/10.1098/rspa.2017.0058

Google Scholar

Crossref

60.

M.

Descroches

,

J.

Guckenheimer

,

B.

Krauskopf

,

C.

Kuehn

,

H.

Osinga

, and

M.

Wechselberger

,

SIAM Rev.

54

,

211

–

288

(

2012

).

https://doi.org/10.1137/100791233

Google Scholar

Crossref

61.

M.

Desroches

,

M.

Krupa

, and

S.

Rodrigues

,

J. Math. Biol.

67

,

989

–

1017

(

2013

).

https://doi.org/10.1007/s00285-012-0576-z

Google Scholar

Crossref

PubMed

62.

R.

Benoît

and

C.

Lobry

,

C. R. Acad. Sci. Paris Ser. I

294

,

483

–

488

(

1982

).

Google Scholar

63.

J.-L.

Callot

,

F.

Diener

, and

M.

Diener

,

C. R. Acad. Sci. Paris Ser. A

286

,

1059

–

1061

(

1978

).

Google Scholar

64.

F.

Diener

and

M.

Diener

,

C. R. Acad. Sci. Paris Ser. I

297

,

577

–

588

(

1983

).

Google Scholar

65.

A.

Roberts

and

P.

Glendinning

,

Chaos

24

,

023138

(

2014

).

https://doi.org/10.1063/1.4885502

Google Scholar

Crossref

PubMed

66.

E. K.

Ersöz

,

M.

Desroches

, and

M.

Krupa

,

Phys. D

349

,

46

–

61

(

2017

).

https://doi.org/10.1016/j.physd.2017.02.016

Google Scholar

Crossref

67.

Y.

Miyagawa

,

A.

Tsujimura

,

K.

Fujita

et al.,

NeuroImage

36

,

830

–

842

(

2007

).

https://doi.org/10.1016/j.neuroimage.2007.03.055

Google Scholar

Crossref

PubMed

68.

J. G.

Bohlen

,

J. P.

Held

, and

M. O.

Sanderson

,

Arch. Sex. Behav.

9

,

503

–

521

(

1980

).

https://doi.org/10.1007/BF01542155

Google Scholar

Crossref

PubMed

69.

C. M.

Meston

and

R. J.

Levin

, “Female orgasm dysfunction,” in Handbook of Sexual Dysfunction, edited by R. Balon and R. T. Segraves (CRC Press, Boca Raton, FL, 2005).

70.

K. R.

Turley

and

D. L.

Rowland

,

BJI Int.

112

,

442

–

452

(

2013

).

https://doi.org/10.1111/bju.12011

Google Scholar

Crossref

71.

L. Q. P.

Paterson

,

E. S.

Jin

,

R.

Amsel

, and

Y. M.

Binik

,

J. Sex Res.

51

,

801

–

813

(

2014

).

https://doi.org/10.1080/00224499.2013.867922

Google Scholar

Crossref

PubMed

72.

T. H. C.

Krüger

,

P.

Haake

,

U.

Hartmann

,

M.

Schedlowski

, and

M. S.

Exton

,

Neurosci. Biobehav. Rev.

26

,

31

–

44

(

2002

).

https://doi.org/10.1016/S0149-7634(01)00036-7

Google Scholar

Crossref

PubMed

73.

T. H. C.

Krüger

,

U.

Hartmann

, and

M.

Schedlowski

,

World J. Urol.

23

,

130

–

138

(

2005

).

https://doi.org/10.1007/s00345-004-0496-7

Google Scholar

Crossref

PubMed

75.

B. R.

Komisaruk

,

B.

Whipple

,

S.

Nasserzadeh

, and

C.

Beyer-Flores

,

The Orgasm Answer Guide

(

John Hopkins University press

,

Baltimore, MD

,

2009

).

Google Scholar

Crossref

76.

W. E.

Bradley

,

G. W.

Timm

,

J. M.

Gallagher

, and

B. K.

Johnson

,

Urology

26

,

4

–

9

(

1985

).

https://doi.org/10.1016/0090-4295(85)90243-2

Google Scholar

Crossref

PubMed

77.

W. O.

Kessler

,

Urol. Clin. North Am.

15

,

81

–

86

(

1988

).

https://doi.org/10.1016/S0094-0143(21)00872-7

Google Scholar

Crossref

PubMed

78.

D. G.

Hatzichristou

,

K.

Hatzimouratidis

,

E.

Ioannides

,

K.

Yannakoyorgos

,

G.

Dimitriadis

, and

A.

Kalinderis

,

J. Urol.

159

,

1921

–

1926

(

1998

).

https://doi.org/10.1016/S0022-5347(01)63197-5

Google Scholar

Crossref

PubMed

79.

T.

Grafke

and

E.

Vanden-Eijnden

,

J. Stat. Mech.

2017

,

093208

.

https://doi.org/10.1088/1742-5468/aa85cb

Crossref

80.

E.

Forgoston

and

I. B.

Schwartz

,

SIAM J. Appl. Dyn. Syst.

8

,

1190

–

1217

(

2009

).

https://doi.org/10.1137/090755710

Google Scholar

Crossref

81.

P. J.

Thomas

and

B.

Lindner

,

Phys. Rev. E

99

,

062221

(

2019

).

https://doi.org/10.1103/PhysRevE.99.062221

Google Scholar

Crossref

PubMed

82.

A.

Pérez-Cervera

,

B.

Lindner

, and

P. J.

Thomas

,

Biol. Cybernet.

116

,

219

–

234

(

2022

).

https://doi.org/10.1007/s00422-022-00929-6

Google Scholar

Crossref

83.

C. W.

Gardiner

,

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

(

Springer

,

Berlin

,

2004

).

Google Scholar

Crossref

84.

P.

Lancaster

and

L.

Rodman

,

Algebraic Riccati Equations

(

Clarendon Press

,

Oxford

,

1995

).

Google Scholar

85.

G. N.

Mil’shtein

and

L. B.

Ryashko

,

J. Appl. Math. Mech.

59

,

47

–

56

(

1995

).

https://doi.org/10.1016/0021-8928(95)00006-B

Google Scholar

Crossref

86.

I.

Bashkirtseva

and

L.

Ryashko

,

Chaos

21

,

047514

(

2011

).

https://doi.org/10.1063/1.3647316

Google Scholar

Crossref

PubMed

87.

I.

Bashkirtseva

and

L.

Ryashko

,

Phys. Rev. E

83

,

061109

(

2011

).

https://doi.org/10.1103/PhysRevE.83.061109

Google Scholar

Crossref

88.

L.

Ryashko

and

I.

Bashkirtseva

,

Math. Model. Nat. Phenom.

10

,

130

–

141

(

2015

).

https://doi.org/10.1051/mmnp/201510209

Google Scholar

Crossref

89.

A. J.

McKane

and

T. J.

Newman

,

Phys. Rev. Lett.

94

,

218102

(

2005

).

https://doi.org/10.1103/PhysRevLett.94.218102

Google Scholar

Crossref

PubMed

90.

D.

Alonso

,

A. J.

McKane

, and

M.

Pascual

,

J. R. Soc. Interface

4

,

575

–

582

(

2007

).

https://doi.org/10.1098/rsif.2006.0192

Google Scholar

Crossref

PubMed

91.

M. I.

Freidlin

and

A. D.

Wentzell

,

Random Perturbations of Dynamical Systems

(

Springer-Verlag

,

Berlin

,

2012

).

Google Scholar

Crossref

92.

T.

Grafke

and

E.

Vanden-Eijnden

,

Chaos

29

,

063118

(

2019

).

https://doi.org/10.1063/1.5084025

Google Scholar

Crossref

PubMed

93.

M.

Heymann

and

E.

Vanden-Eijnden

,

Phys. Rev. Lett.

100

,

140601

(

2007

).

https://doi.org/10.1103/PhysRevLett.100.140601

Google Scholar

Crossref

94.

M.

Heymann

and

E.

Vanden-Eijnden

,

Commun. Pure Appl. Math.

61

,

1052

–

1117

(

2008

).

https://doi.org/10.1002/cpa.20238

Google Scholar

Crossref

95.

A. J.

Bray

and

A. J.

McKane

,

Phys. Rev. Lett.

62

,

493

–

496

(

1989

).

https://doi.org/10.1103/PhysRevLett.62.493

Google Scholar

Crossref

PubMed

96.

S.

Beri

,

R.

Mannella

,

D. G.

Luchinsky

,

A. N.

Silchenko

, and

P. V. E.

McClintock

,

Phys. Rev. E

72

,

036131

(

2005

).

https://doi.org/10.1103/PhysRevE.72.036131

Google Scholar

Crossref

97.

Z.

Xiang

,

W.

Ren

, and

W.

E

,

J. Chem. Phys.

128

,

104111

(

2008

).

Google Scholar

PubMed

98.

Z.

Xiang

and

W.

E

,

Commun. Math. Sci.

8

,

341

–

355

(

2010

).

Google Scholar

99.

Q.

Yu

and

X.

Liu

,

Chaos

31

,

093110

(

2021

).

https://doi.org/10.1063/5.0056784

Google Scholar

Crossref

PubMed

100.

T.

Grafke

,

T.

Schäfer

, and

E.

Vanden-Eijnden

, “Long term effects of small random perturbations on dynamical systems: Theoretical and computational tools,” in Recent Progress and Modern Challenges in Applied Mathematics, Modeling and Computational Science, Fields Institute Communications, edited by R. Melnik, R. Makarov, and J. Belair (Springer, New York, 2017), pp. 17–55.

101.

D.

Dahiya

and

M.

Cameron

,

J. Sci. Comp.

75

,

1351

–

1384

(

2018

).

https://doi.org/10.1007/s10915-017-0590-9

Google Scholar

Crossref

102.

M. K.

Cameron

,

Phys. D

241

,

1532

–

1550

(

2012

).

https://doi.org/10.1016/j.physd.2012.06.005

Google Scholar

Crossref

103.

D.

Dahiya

and

M.

Cameron

,

Phys. D

382–383

,

33

–

45

(

2018

).

https://doi.org/10.1016/j.physd.2018.07.002

Google Scholar

Crossref

104.

Y.

Li

,

J.

Wang

, and

X.

Liu

,

Commun. Nonlinear Sci. Numer. Simul.

91

,

105441

(

2020

).

https://doi.org/10.1016/j.cnsns.2020.105441

Google Scholar

Crossref

105.

Q.

Yu

and

X.

Liu

,

J. Stat. Mech.

2022

,

013207

.

https://doi.org/10.1088/1742-5468/ac3e73

Crossref

106.

Y.

Li

, see https://github.com/liyangnuaa/Quasi-threshold-phenomenon-in-noise-driven-Higgins-model for “Quasi-Threshold Phenomenon in Noise-Driven Higgins Model.”

107.

R. S.

Maier

and

D. L.

Stein

,

SIAM J. Appl. Math.

57

,

752

–

790

(

1997

).

https://doi.org/10.1137/S0036139994271753

Google Scholar

Crossref

108.

M.

Sipski

,

C. J.

Alexander

, and

O.

Gómez-Marín

,

Spinal Cord

44

,

798

–

804

(

2006

).

https://doi.org/10.1038/sj.sc.3101954

Google Scholar

Crossref

PubMed

109.

M.

Alexander

and

R. C.

Rosen

,

J. Sex Marital Ther.

34

,

308

–

324

(

2008

).

https://doi.org/10.1080/00926230802096341

Google Scholar

Crossref

PubMed

110.

S. M.

Castle

,

L. C.

Jenkins

,

E.

Ibrahim

,

T. C.

Aballa

,

C. M.

Lynne

, and

N. L.

Brackett

,

Spinal Cord

52

,

S27

–

S29

(

2014

).

https://doi.org/10.1038/sc.2014.110

Google Scholar

Crossref

PubMed

111.

M.

Sipski

,

C. J.

Alexander

,

O.

Gómez-Marín

, and

J.

Spalding

,

J. Urol.

177

,

247

–

251

(

2007

).

https://doi.org/10.1016/j.juro.2006.08.059

Google Scholar

Crossref

PubMed

112.

W. G.

Dail

,

G.

Walton

, and

M. P.

Olmsted

,

J. Autonom. Nerv. Syst.

28

,

251

(

1989

).

https://doi.org/10.1016/0165-1838(89)90153-7

Google Scholar

Crossref

113.

B. D.

Sachs

,

Neurosci. Biobehav. Rev.

19

,

211

(

1995

).

https://doi.org/10.1016/0149-7634(94)00063-7

Google Scholar

Crossref

PubMed

114.

B. S.

Matthews

and

J. W.

Wells

,

J. Sex Educ. Ther.

9

,

26

–

31

(

1983

).

https://doi.org/10.1080/01614576.1983.11074778

Google Scholar

Crossref

115.

Y.

Hovav

,

M.

Dan-Goor

,

H.

Yaffe

, and

M.

Almagor

,

Fertil. Steril.

72

,

364

–

365

(

1999

).

https://doi.org/10.1016/S0015-0282(99)00239-3

Google Scholar

Crossref

PubMed

116.

A.

Gul

,

E.

Yuruk

, and

E. C.

Serefoglu

,

Rev. Int. Androl.

18

,

21

–

26

(

2020

); available at https://www.sciencedirect.com/science/article/abs/pii/S1698031X18300815.

Google Scholar

PubMed

117.

C. K.-C.

Yu

and

W.

Fu

,

Dreaming

21

,

197

–

212

(

2011

).

https://doi.org/10.1037/a0024085

Google Scholar

Crossref

118.

J. M.

Davidson

, “The psychobiology of sexual experience,” in The Psychobiology of Consciousness, edited by J. M. Davidson and R. J. Davidson (Plenum Press, New York, 1980), pp. 271–332.

119.

A.

Safron

,

Socioaffect. Neurosci. Psychol.

6

,

31763

(

2016

).

https://doi.org/10.3402/snp.v6.31763

Google Scholar

Crossref

PubMed

120.

B. R.

Komisaruk

, and

B.

Whipple

,

Sex. Relatsh. Ther.

26

,

356

–

372

(

2011

).

https://doi.org/10.1080/14681994.2011.649252

Google Scholar

Crossref

121.

D.

Sylva

,

A.

Safron

,

A. M.

Rosenthal

,

P. J.

Reber

,

T. B.

Parrish

, and

J. M.

Bailey

,

Horm. Behav.

64

,

673

–

684

(

2013

).

https://doi.org/10.1016/j.yhbeh.2013.08.003

Google Scholar

Crossref

PubMed

122.

A.

Koseska

,

E.

Volkov

, and

J.

Kurths

,

Phys. Rep.

531

,

173

–

199

(

2013

).

https://doi.org/10.1016/j.physrep.2013.06.001

Google Scholar

Crossref

123.

A.

Pikovsky

,

M.

Rosenblum

, and

J.

Kurths

,

Synchronization: A Universal Concept in Nonlinear Sciences

(

Cambridge University Press

,

Cambridge

,

2001

).

Google Scholar

Crossref

124.

S.

Majhi

,

B. K.

Bera

,

D.

Ghosh

, and

M.

Perc

,

Phys. Life Rev.

28

,

100

–

121

(

2019

).

https://doi.org/10.1016/j.plrev.2018.09.003

Google Scholar

Crossref

PubMed

125.

M.

Apter

,

The Dangerous Edge: The Psychology of Excitement

(

The Free Press

,

New York

,

1992

).

Google Scholar

126.

M. J.

Apter

,

Zigzag: Reversal and Paradox in Human Personality

(

Troubador Publishing Ltd.

,

Leicester

,

2018

).

Google Scholar

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

I. INTRODUCTION

II. AN ORGASM IS NOT A CATASTROPHE

III. MODEL DERIVATION

IV. STABILITY ANALYSIS AND NUMERICAL SIMULATIONS

V. STOCHASTIC EFFECTS

A. Stochastic oscillations near equilibrium

B. Large deviations: An orgasm as an escape

VI. DISCUSSION

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

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

I. INTRODUCTION

II. AN ORGASM IS NOT A CATASTROPHE

III. MODEL DERIVATION

IV. STABILITY ANALYSIS AND NUMERICAL SIMULATIONS

V. STOCHASTIC EFFECTS

A. Stochastic oscillations near equilibrium

B. Large deviations: An orgasm as an escape

VI. DISCUSSION

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only