Multilevel Mathematical Model for Manual Therapy: Integration of Viscoelastic Mechanics, Neurophysiology of Pain and Clinical Dose-response

Methods

General design of the model

The model is organized into three cascaded levels (Figure 1): the output of Level 1 (afferent firing rate) becomes the input of Level 2 (neurophysiological dynamics), whose output feeds into the integration function of Level 3 (clinical dose-response). This architecture reflects the actual physiological causal chain and allows each level to be calibrated or replaced independently, a desirable property in models intended to evolve with empirical data.

Figure 1: Simulated trajectories of the three levels of the model during a 60-second Manual Therapy session (shaded region). From top to bottom: applied strain ε(t), tissue stress σ(t), afferent rate R_aff(t), neurophysiological dynamics T(t)/S(t), effective dose D_ef(t), and predicted clinical response ΔPain(t).

Level 1 — Tissue Mechanics and Mechanotransduction

Standard Linear Solid Model: Soft tissue is modeled as a standard linear solid (SLS), equivalent to the Zener model (spring in parallel with a Maxwell element), whose constitutive equation is:

σ(t) + τσ · σ̇(t) = E₀ · ε(t) + E₀τε · ε̇(t)                          (Eq. 1)

where σ(t) is the stress [kPa], ε(t) the strain, E₀ the instantaneous elastic modulus, and τσ, τε the relaxation and retardation time constants, respectively. The SLS was selected over single-element models (Kelvin-Voigt or Maxwell) for its ability to simultaneously reproduce stress relaxation and retarded strain, behaviors documented in musculoskeletal soft tissue [12,13].

For a sustained slip technique (step deformation ε₀), the analytical solution is:

σ(t) = ε₀ [E∞ + (E₀ − E∞) · e^(−t/τσ)]                       (Eq. 2)

The E∞/E₀ ratio is clinically interpretable as an index of tissue relaxability, with lower values in fibrotic tissues or those with chronic mobility restrictions.

Mechanotransduction: rapidly and slowly adapting populations: Tissue deformation activates cutaneous and subcutaneous mechanoreceptors. Consistent with the functional classification of primary afferents described by Niimi et al. [14] and the evidence on low-threshold mechanoreceptor properties documented in the sensory neurophysiology literature [15], the model distinguishes two populations:

a) Rapidly adapting receptors (analogous to Pacini/Meissner corpuscles): sensitive to the strain rate |ε̇|, predominant in oscillatory techniques. b) Slowly adapting receptors (analogous to Ruffini corpuscles/Merkel discs): sensitive to the magnitude of sustained strain |ε|, predominant in sustained slip.

Both populations are modeled with a saturable Hill function (Eq. 3-4), and the net afferent signal is its point maximum (Eq. 5):

R_fast(t) = R_max · | ε ̇(t)|ⁿ / (K_m,r ⁿ + | ε ̇(t)|ⁿ)        (Eq. 3)

R_slow(t) = R_max · | ε (t)|ⁿ / (K_m,l ⁿ + | ε (t)|ⁿ)     (Eq. 4)

R_aff(t) = max (R_fast(t), R_slow(t))                                 (Eq. 5)

Level 2 — Neurophysiology: spinal gating and central sensitization

Transmission cell system: The Gate Control Theory of pain [5] is the neurophysiological framework with the strongest experimental support for the analgesic effects of TM [1,6]. Following the mathematical formalization of Britton and Skevington [16] and the computational extension of the dorsal horn developed by Crodelle and Maia [9], the activity of the transmitting cell (T cell) of the dorsal horn is described by the ODE:

dT/dt = −α · T + β · R_noci(t) − γ · R_aff(t) · G(T) + δ · D(t)                                (Eq. 6)

where T(t) is the correlate of the ascending nociceptive message, R_noci(t) the basal nociceptive afferent (Aδ/C fibers, independent of the technique), G(T) the presynaptic inhibition of the substantia gelatinosa, D(t) the descending cortical control (expectation, therapeutic alliance [6] ), and α the natural decay rate.

The function G(T) is modeled as a sigmoid (Eq. 7), reproducing the non-linearity of presynaptic inhibition:

G(T) = 1 / (1 + e^(−k · (T − T₀)))                                 (Eq. 7)

This formulation is mathematically equivalent to that interpreted by Bhatt et al. [10] in their multiscale dorsal horn model, which also incorporates conductance-based units and produces firing patterns consistent with experimental data.

Central sensitization dynamics: Central sensitization—amplification of the nociceptive signal by synaptic plasticity in the dorsal horn— [7,10] is modeled as a second state variable S(t) that evolves on a much slower timescale than T(t):

dS/dt = ε · [T(t) − S(t)], ε ≪ α                                       (Eq. 8)

The temporal hierarchy ε ≪ α produces, in an emergent manner —without the need for ad hoc conditionals— the lower speed of therapeutic response observed in chronic pain versus acute pain, a well-documented clinical pattern that until now lacked quantitative expression in the context of TM.

Level 3 — Clinical Dose-Response

The effective dose D_ef(t) is defined as the cumulative time average of the net mechanoreceptive balance (Eq. 9), normalized by time to make it independent of the gross session duration:

D_ef(t) = (1/t) · ∫₀ᵗ [w₁ · R_aff(τ) − w₂ · T(τ)]₊ dτ    (Eq. 9)

The clinical change in pain is predicted by a sigmoidal Emax curve (Eq. 10), standard in clinical pharmacology and transferable to the concept of 'mechanical dose':

ΔPain(t) = Emax · D_ef(t)ʰ / (ED₅₀ʰ + D_ef(t)ʰ)     (Eq. 10)

All three parameters have direct clinical interpretation: Emax is the therapeutic ceiling (dependent on chronicity and comorbidity), ED₅₀ is the patient's sensitivity to the intervention, and h determines whether the response is gradual (h<2) or an abrupt threshold (h>3). This parameterization is compatible with the recommendation of the dose-response models in physiotherapy reviewed by Bhogal et al. [17].

Numerical implementation and validation protocol

The system of ODEs (Eqs. 6-8) was solved numerically using the adaptive Runge-Kutta integrator RK45 (SciPy, Python 3.11), with tolerances rtol=10⁻⁶ and atol=10⁻⁸ to ensure convergence. The floor function of the system (T ≥ 0) was implemented using a smooth logistic sigmoid G_floor(T) = 1/(1+e^(−T)) that preserves the C¹ continuity necessary for adaptive integration, avoiding discontinuities that would generate infinitesimal steps.

The physiological coherence validation protocol included: (i) sensitivity analysis to the intensity of the technique (ε₀ ∈ [0.005–0.18]); (ii) sensitivity analysis to duration; (iii) comparison between sustained step vs. oscillatory techniques (2 Hz, similar to Maitland grade III [18n19]); (iv) effect of the patient's baseline irritability on the response; and (v) analysis of dose-response curves for three chronicity profiles (acute, subacute, complex chronic).

Justification of the parameterization and structural sensitivity

The numerical values assigned to the constants in Equations 6–8 (α, β, γ, δ, k, T 0 , ε) and to the Level 1 parameters (E 0 , E∞, τσ, τε, K m,r , K m,l , n) were selected by order of magnitude from the physiological ranges reported in the cited literature [12,14,15] and not by fitting to patient data. The separation of time scales ε ≪ α (Eq. 8) is not an arbitrary assumption: it is the minimum condition necessary to reproduce, without additional ad hoc parameters, the difference of an order of magnitude observed between the fast kinetics of spinal gating (T, τ≈1.25 s) and the slow kinetics of central sensitization (S, τ≈50 s) clinically documented [7].

To assess the robustness of the results described in sections 3.1–3.3 against the uncertainty in this parameterization, an additional univariate sensitivity analysis was performed, varying each constant by ±20% from its reference value. The percentage closure of the spinal gate (85.3% of the baseline value) varied between 78% and 91%, and the time constant of S(t) varied proportionally with 1/ε, without altering the qualitative ordering of the results in Table 1. This suggests that the qualitative conclusions of the model are reasonably robust to the uncertainty in the current parameterization, although, as discussed in section 4.3, their quantitative accuracy depends on future calibration with real clinical data.

Results

Dynamics of the three levels

Simulation of a standard session (60 s of oscillatory technique at 2 Hz, ε₀=0.08, initiated at t=10 s) produced physiologically coherent trajectories at all levels (Figure 1).

The oscillatory technique generated bursts of mechanoreceptive afferent input with a maximum R_aff of 98.8 Hz, consistent with firing rate ranges reported in low-threshold mechanoreceptors under mechanical stimulation [14,15].

T-cell activity decreased from its baseline value (T₀ ≈ 11.3) to values close to zero during the application of the technique, with an 85.3% reduction from baseline—an effective closure of the spinal gate—without requiring special adjustment parameters. Upon removal of the stimulus (t=70 s), T rebounded to a value slightly above baseline and then decayed exponentially toward the new equilibrium state, consistent with the post-inhibition rebound mechanism documented in inhibitory interneurons of the dorsal horn [10].

The central sensitization S(t) showed a markedly slower dynamic (τ ≈ 50 s vs τ ≈ 1.25 s for T), with a difference of an order of magnitude between time constants, validating the hierarchical time structure of the model.

Parametric sensitivity analysis

Table 1 summarizes the results of the sensitivity analysis for the dosage parameters.

The intensity-response relationship was monotonic and showed clear saturation above ε₀ = 0.12, consistent with the saturation of the Hill function of the mechanoreceptors. The step/oscillatory comparison showed a difference of 9.2% in D_ef and 8.8% in ΔPain, with a slight advantage for the oscillatory technique, consistent with evidence that positions Maitland grade III-IV techniques as superior in pain reduction under conditions of high irritability [18,19].

The effect of baseline irritability was statistically expected but quantitatively moderate (−11.3% in ΔPain when quadrupling R_noci), consistent with the clinical observation that TM maintains efficacy even in cases of high irritability when appropriate techniques and doses are selected.

Dose-response curves according to chronicity

Figure 2 shows the differential behavior of the parameters Emax, ED₅₀, and h in the dose-response curves for the three chronicity profiles. The three parameters were independently identifiable (orthogonal gradients in response space), confirming that the model is parametrically identifiable—a necessary prerequisite for estimating parameters from longitudinal clinical data.

Figure 2: Dose-response curves: Emax. Left: effect of the ED₅₀ parameter (patient sensitivity). Center: effect of the Hill coefficient h (graduality of response). Right: effect of Emax (therapeutic ceiling) for three chronicity profiles.

Corrections identified during validation

The validation process identified and corrected two design errors that would have compromised the physiological consistency of the model:

(a) Artificial saturation of the effective dose: the initial formulation of D_ef(t) as a cumulative integral without normalization produced values that saturated the Emax curve across virtually the entire physiological range, eliminating parametric sensitivity. The correction (temporal normalization, Eq. 9) restored the graded dose-response relationship.

(b) Insensitivity to sustained gliding: the initial formulation of R_aff based exclusively on |ε̇| predicted an almost zero response for step-like techniques, contradicting the clinical evidence of the efficacy of sustained gliding [19,20]. The incorporation of a second population of slowly adapting mechanoreceptors (Eq. 4), sensitive to |ε|, aligned the model with known sensory physiology.

Practical guide for selecting a technique according to clinical objective

The results of sections 3.2 and 3.3 allow us to formulate a preliminary guide, subject to clinical validation, on which model configuration (type of technique, dose and chronicity profile) is most appropriate according to the therapeutic objective sought (Table 2).

Discussion

The proposed model represents the first documented attempt to integrate, within a computable system of coupled equations, the three levels of action of Manual Therapy: tissue mechanics, the neurophysiology of pain, and the clinical dose-response relationship. Its main contribution is not to accurately predict the response of a specific patient—which requires calibration with real data—but rather to offer a formal architecture that makes explicit the functional relationships between dosage variables and outcome variables, which have been implicit until now in the qualitative conceptual frameworks of Manual Therapy.

Theoretical support of the three levels

Level 1 (tissue mechanics) employs the SLS, whose superiority over single-element models for describing the viscoelastic behavior of soft tissues under constant and transient loads was confirmed by Ahmed et al. [11] in corneal tissue, and by Holzapfel et al. [12] in the general constitutive framework of soft tissues. The extension to two mechanoreceptor populations (rapid/slow adaptation) is supported by the electrophysiological characterization of Niimi et al. [14], which experimentally distinguishes between Nav1.8-positive (nociceptive) and Nav1.8-negative (low-threshold mechanoreceptor) subtypes by their adaptation pattern—SA or RA response—and their mechanical thresholds.

Level 2 (neurophysiology) connects directly with the line of research on computational modeling of the dorsal horn formally initiated by Britton and Skevington [16], extended by Bhatt et al. [10] with multiscale conductance models, and framed within the systematic review by Lang et al. [8] that identified mechanistic models of nociceptive processing as the category with the greatest potential for generating testable clinical hypotheses. The dynamics of central sensitization (Eq. 8) are consistent with the spinal LTP and NMDA-dependent mechanisms described by Woolf [7] and with the timescale differences between acute and chronic pain observed clinically.

Level 3 (dose-response) applies the sigmoidal Emax structure, standard in clinical pharmacology (Sheiner et al. model), to the 'mechanical dose' domain. This translation is analogous to that used in the dose-response review of spasticity in physiotherapy [17] and has been implicitly assumed—without formalization—in the practice of grading Maitland techniques (grades I–IV) as a dose proxy.

Implications for research in Manual Therapy

The model generates testable quantitative hypotheses that are not deducible from current qualitative frameworks. For example: (i) the difference in response between oscillatory and sustained techniques should be approximately 9-12% in acute patients with normal tissue mobility, but could increase in fibrotic tissues (high E∞/E₀) where slow adaptation is attenuated; (ii) in chronic pain with high central sensitization (high S), the response to a single session should be markedly lower than in acute pain—not proportional—because D_ef is attenuated by the higher baseline level of T; (iii) the descending control D(t) should act multiplicatively on the response (not additively), predicting specific interactions between therapeutic alliance and technique intensity.

These hypotheses can be directly designed as observational or experimental studies with standard outcome variables (EVA, PPT, ROM), closing the loop between theoretical modeling and empirical evidence.

Limitations

Four main limitations should be considered: (i) The model parameters have not been calibrated with real patient data; the values used are reasonable by order of magnitude and consistent with physiological ranges, but their current quantitative precision is limited. At this stage, a formal estimate of the uncertainty (confidence intervals or standard error) associated with these values is not available, since they are derived from literature ranges and not from a statistical fit to a sample of patients; it is recommended that the future calibration phase explicitly report this parametric uncertainty, for example, through bootstrap resampling or Bayesian analysis on the collected clinical data. (ii) The top-down control D(t) reduces a complex psychosocial variable (expectation, context, therapeutic alliance, catastrophizing) to a scalar, sacrificing relevant nuances. (iii) The model does not incorporate cumulative effects between sessions or local inflammatory processes. (iv) The validation performed is for internal consistency, not predictive accuracy. The calibration and clinical validation roadmap should follow the steps proposed in the research agenda on TM mechanisms by Bialosky et al.

Conclusions

1. The first multilevel mathematical model is proposed that connects, in a computable way and with clinically interpretable parameters, the force applied during a Manual Therapy technique with the predicted change in pain.
2. The three-level coupled architecture (viscoelastic SLS → spinal gate EDOs → Emax curve) produces physiologically consistent behavior and is parametrically identifiable, conditions necessary for its calibration with real clinical data.
3. The distinction between rapidly and slowly adapting mechanoreceptor populations (Hill differentiated on |ε̇| and |ε|, respectively) is necessary for the model to predict substantial responses for both oscillatory techniques and sustained glides.
4. The central sensitization variable S(t) with a slow time scale (ε ≪ α) emergently produces the lowest therapeutic response speed in chronic pain documented clinically, without requiring ad hoc parameterization.
5. In its current state, the model provides a formal framework for generating testable hypotheses and designing clinical support tools, but it is not a validated clinical prediction tool. Its direct clinical application is contingent upon parametric calibration with real patient data, followed by prospective validation.

Data and code availability

The complete code for the numerical implementation (Python 3.11, SciPy, RK45 integrator with rtol=10⁻⁶ and atol=10⁻⁸), the parametric values used in each simulation, and the scripts used to generate Figure 1, Figure 2, and Table 1 will be deposited in an open-access public repository (e.g., Zenodo or OSF) at the time of publication, instead of being available upon request to the author. This will allow for the complete reproduction of the simulations and sensitivity analysis described in sections 2.5, 2.6, and 3.2 without requiring direct contact with the respective author.