Nonlinear dynamics in pulse-modulated feedback drug dosing

TL;DR

This work analyzes pulse-modulated feedback for dosing a neuromuscular blocker, showing that a hybrid closed-loop—combining a linear PK block with a nonlinear Hill PD and impulsive control—can exhibit rich nonlinear dynamics. The authors formulate a discrete map $X_{n+1}=Q(X_n)$ to capture the inter-firing dynamics, and perform bifurcation analysis to reveal 1-cycle, multistable, and chaotic regimes. A design procedure selects piecewise-affine modulation functions $\Phi$ and $F$ to realize a robust 1-cycle under nominal PKPD parameters, while simulations reveal how parameter variations can lead to 2-cycle or 3-cycle attractors and, in some regions, deterministic chaos. The study emphasizes that safe, effective closed-loop dosing requires bifurcation-informed design to avoid undesirable attractors and ensure patient safety. All mathematical notation is expressed with explicit $...$ delimiters to preserve interpretability in downstream analyses.

Abstract

Pulse-modulated feedback is utilized in drug dosing to mimic sustained over a longer period of time manual discrete dose administration, the latter is in contrast with continuous drug infusion. The intermittent mode of dosing calls for a hybrid (continuous-discrete) modeling of the closed-loop system, where the pharmacokinetics and pharmacodynamics of the drug are captured by differential equations whereas the control law is described by difference equations. Hybrid dynamics are highly nonlinear which complicates formal design of pulse-modulated feedback. This paper demonstrates complex nonlinear dynamical phenomena arising in a simple control system of dosing a neuromuscular blockade agent in anesthesia. Along with the nominal periodic regimen, undesirable nonlinear behaviors, i.e. periodic solutions of high multiplicity, multistability, as well as deterministic chaos, are shown to exist. It is concluded that design of feedback drug dosing algorithms based on a hybrid paradigm has to be informed by a thorough bifurcation analysis in order to secure patient safety.

Figures (9)

• Figure 1: Convergence to the 1-cycle from $x(0)=0$ in the NMB model with $(\bar{\alpha},\bar{\gamma})$ stabilized by the modulation function slopes $F^\prime(y_0)=-0.15$, $\Phi^\prime(y_0)=0.29$. Top plot: the nonlinear output $y(t)$. The horizontal black dashed lines mark $\inf_t y(t)$ and $\sup_{t\in \lbrack T,5T\rbrack} y(t)$. The stationary output corridor values for the 1-cycle are marked in red. Bottom plot: the linear output $\bar{y}(t)$.
• Figure 2: Convergence to the 2-cycle from $x(0)=0$ in the NMB model with $(\alpha=0.11, \bar{\gamma})$ stabilized by the modulation function slopes $F^\prime(y_0)=-0.15$, $\Phi^\prime(y_0)=0.29$. Top plot: the nonlinear output $y(t)$. Bottom plot: the linear output $\bar{y}(t)$. The desired output corridor values for the 1-cycle are marked in red. The horizontal black dashed lines mark $\inf_t y(t)$ and $\sup_{t\in \lbrack T,5T\rbrack} y(t)$.
• Figure 3: Bifurcation diagram for $F_1=150.0$, $F_2=400.0$, $\Phi_1=11.0$, $\Phi_2=50.0$. Here $\varphi(\bar{y}_0)=2.1256$ for $\bar{y}_0=13.6249$ and $\alpha=\alpha_*$, $\alpha_*\approx0.0374\approx \bar{\alpha}$. The intervals $0.02<\alpha<\alpha_\mathrm{Flip}$ and $\alpha\gtrapprox\alpha_\mathrm{Flip}^\mathrm{BCB}$ are the regions of the stability of a fixed point. Here $\alpha_\mathrm{Flip}\approx 0.107$, $\alpha_\mathrm{Flip}^\mathrm{BCB}\approx 0.277$. The dotted line marks the unstable fixed point. Between the points $\alpha_\mathrm{Flip}$ and $\alpha_\mathrm{Flip}^\mathrm{BCB}$ exists the stable 2-cycle.
• Figure 4: 2D bifurcation diagram for $F'(y_0)=-1.0$, and $\Phi'(y_0)=4.0$, $F_1=150.0$, $F_2=400.0$, $\Phi_1=5.0$, $\Phi_2=50.0$, $0.027<\alpha<0.0524$, $1.403<\gamma<5.5619$, $a_1=-0.0374$, $a_2=-0.1496$, $a_3=-0.374$, $\mathrm{g}_1=0.0374$, $\mathrm{g}_2=0.056$, $T=20.0$, $\lambda=300.0$.
• Figure 5: 1D bifurcation diagram for $F'(y_0)=-1.0$, and $\Phi'(y_0)=4.0$, $F_1=150.0$, $F_2=400.0$, $\Phi_1=5.0$, $\Phi_2=50.0$, $\gamma=2.6677$, $0.027<\alpha<0.0524$, $a_1=-0.0374$, $a_2=-0.1496$, $a_3=-0.374$, $\mathrm{g}_1=0.0374$, $\mathrm{g}_2=0.056$, $T=20.0$, $\lambda=300.0$, $k_1=40.87311$, $k_2=-9.819844$, $k_3=294.7817$, $k_4=2.45496$, $\varphi(\bar{y}_0)=2.125605$, $\varphi'(\bar{y}_0)=-0.4073$. The multistability region is $\alpha >\alpha_\mathrm{SN}$ (in yellow).