Hybrid statistical and mechanistic mathematical model guides mobile health intervention for chronic pain

Hybrid statistical and mechanistic mathematical model guides mobile health intervention for chronic pain

Sara M. Clifton sclifton@u.northwestern.edu Department of Engineering Sciences and Applied Mathematics, Northwestern University, Evanston, IL 60208, USA    Chaeryon Kang Department of Biostatistics, University of Pittsburgh, Pittsburgh, PA 15261, USA    Jingyi Jessica Li Department of Statistics, University of California Los Angeles, Los Angeles, CA 90095, USA    Qi Long Department of Biostatistics and Epidemiology, University of Pennsylvania, Philadelphia, PA 19104, USA    Nirmish Shah Department of Medicine, Duke University, Durham, NC 27710, USA    Daniel M. Abrams Department of Engineering Sciences and Applied Mathematics, Northwestern University, Evanston, IL 60208, USA Northwestern Institute for Complex Systems, Northwestern University, Evanston, IL 60208, USA Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA
July 5, 2019

Nearly a quarter of visits to the Emergency Department are for conditions that could have been managed via outpatient treatment; improvements that allow patients to quickly recognize and receive appropriate treatment are crucial. The growing popularity of mobile technology creates new opportunities for real-time adaptive medical intervention, and the simultaneous growth of “big data” sources allows for preparation of personalized recommendations. Here we focus on the reduction of chronic suffering in the sickle cell disease community. Sickle cell disease is a chronic blood disorder in which pain is the most frequent complication. There currently is no standard algorithm or analytical method for real-time adaptive treatment recommendations for pain. Furthermore, current state-of-the-art methods have difficulty in handling continuous-time decision optimization using big data. Facing these challenges, in this study we aim to develop new mathematical tools for incorporating mobile technology into personalized treatment plans for pain. We present a new hybrid model for the dynamics of subjective pain that consists of a dynamical systems approach using differential equations to predict future pain levels, as well as a statistical approach tying system parameters to patient data (both personal characteristics and medication response history). Pilot testing of our approach suggests that it has significant potential to predict pain dynamics given patients’ reported pain levels and medication usages. With more abundant data, our hybrid approach should allow physicians to make personalized, data driven recommendations for treating chronic pain.

I Introduction

In the fields of physics, chemistry, and engineering, models are often derived from mechanistic fundamental laws expressed in the form of differential equations. Resulting “dynamical systems” models can be used both to gain intuition into the expected behavior of the system, and to make specific predictions about results of experiments (e.g., see strogatz2014nonlinear ()). In fields such as social sciences, bioinformatics, and medicine, models are often constructed from data via statistical inference, without direct derivation from fundamental principles (e.g. freedman2009statistical ()). The mechanistic and statistical approaches to mathematical modeling have different advantages. The former allows prior knowledge to be introduced and validated or rejected based on the success of the model. The latter requires almost no a-priori information about how the system is expected to behave.

In this paper, we present a hybrid approach to mathematical modeling that incorporates both mechanistic and statistical elements, with the goal of gaining a deeper understanding of the human experience of subjective pain. Specifically, we hope to predict how patient-reported pain levels vary over time based on medication dosage information and other patient characteristics.

Figure 1: Smartphone app. Sample images of SMART app for iPhone/Android smartphone devices.

i.1 Application to pain

Sickle cell disease (SCD) is a chronic illness associated with frequent medical complications and hospitalizations. Approximately 90% of acute care visits are for pain events, and 30-day hospital re-utilization rates are alarmingly high platt1991pain (). While factors influencing these high re-utilization rates are poorly understood, close follow-up and continued use of pain medications has been shown to decrease re-hospitalization rates. Mobile technology has become an integral part of health care management, and our recently self-developed mobile application (Sickle cell Mobile Application to Record symptoms via Technology, or SMART app—see Figure 1) for SCD assists with documentation and intervention of pain.

Pain in particular is difficult to quantify and has never before been monitored at the temporal scale we report here across so many patients. It is known that subjective pain, though indeed subjective, is correlated with objective measurable stimuli qualities in experiments (see, e.g., granovsky2008objective (); hughes2002assessment (); stevens1958scale ()). Thus there is reason to believe that subjective pain may follow understandable dynamics in time, especially when mitigated by opioid or non-opioid drugs. Our approach to the problem is motivated by the hope that a reasonable model for pain dynamics will yield some level of predictive power, despite the clear expectation that there will also be significant noise within and across patients. We can attribute the stochastic variation to sources like patient mood, temporal changes in patient state, weather, etc. In contrast, we hope that patient attributes like age, gender, SCD disease type, etc. will remain roughly constant on the time scale of the experiment and allow us to explore possible correlation of these attributes with model parameters.

i.2 Data source: mobile health application

We seek to understand the temporal dynamics of chronic pain as experienced by SCD patients. To that end, we have developed a mobile phone application that allows patients to record medication usage and subjective pain levels (measured on a 0-10 scale) in real time jonassaint2015usability (); shah2014patients ().

Figure 1 shows several images of the app interface, while Figure 2 shows a typical data set resulting from a single patient’s use of the app over the course of several weeks.

Figure 2: Sample pain and medication data from a single patient. Upper panel: patient reported pain (black circles) and model fit (red solid line); red shading indicates model fit plus/minus one standard deviation. Lower panel: long-acting methadone (red solid line) and short-acting oxycodone (blue dashed line) medication concentrations in patient bloodstream as inferred from medication usage reported via the SMART application.

Ii Materials and Methods

ii.1 Data

As of October 2016, data were available from 47 patients using the SMART app. Data sets from 8 of those patients were excluded because of excessive sparsity based on the following criteria:(1) total number of reports ; or (2) pain reports never exceeded zero during the period under consideration. See Table 1 for demographic details of included patients. We denote the sample size .

Demographic characteristics
                                                           A 14
                                                           B 17
                                                           C 8
                                                           Male 16
                                                           Female 23
Age at baseline (years)
                                                           18-34 24
SCD disease type
                                                          Hemoglobin SC 8
                                                          Hemoglobin SS 22
                                                          Hemoglobin SB (Beta) Thalassemia 5
                                                          Beta-Zero Thalassemia 3
                                                          SOAra 1
Hydroxyurea user 27
Folic acid vitamin user 26
Long-acting opioid user 29
Short-acting opioid user 35
Non-opioid user 29

Mean SD (Min, Max)
Number of pain reports 67.2 ( 9.0, 257.0 )
Days of pain reports 164.6 ( 10.3, 435.1 )
Within-patient average VAS score 4.7 ( 0.3, 9.4 )

Mean SD (Min, Max)
Number of pain reports (first 2 weeks) 13.2 ( 2.0, 45.0 )
Number of long-acting opioid doses (first 2 weeks) 6.0 ( 0.0, 35.0 )
Number of short-acting opioid doses (first 2 weeks) 7.2 ( 0.0, 35.0 )
Number of non-opioid doses (first 2 weeks) 2.1 ( 0.0, 12.0 )

Table 1: Patient demographic information and the number of pain reports supplied by patients across entire study.

ii.2 Predictive model

In order to develop a hybrid model that incorporates both a mechanistic a-priori knowledge-driven component and a statistical data-driven component, we divide tasks into two disjoint sets that fit these two categories; see the Discussion section for more context. We begin with a “dynamical systems” model for subjective pain motivated by the hypothesis that human sensory systems function on a roughly “return to setpoint” basis mcruer1974mathematical (); fors1988ability (); britton1995mathematical (); stepan2009delay (). Any model of human pain response, however, will inevitably require specification of a variety of parameters determining the time scale(s) and degree of severity of the response. The statistical modeling tasks employ patient data to infer parameters (1) from patient characteristics and population distributions and (2) from patient-specific pain and medication response history.

To make this more concrete, in Figure 3 we present a flow chart summarizing our approach to the hybrid modeling problem. Steps and A comprise the statistical modeling component; steps B and C comprise the mechanistic modeling component. A further optimization step D builds on the predictions of the hybrid model to allow for a balance between competing demands of pain reduction and medication usage minimization. This paper details steps , , A and B. We leave the remaining steps for future work.

Figure 3: Schematic flowchart showing model framework. Rounded rectangles represent modeling or computation steps, rhombuses represent data inputs or outputs, and diamond represents decision step. Items and are only necessary for initialization of the model. Items A through E are the focus of this paper.

ii.2.1 Mechanistic component (for every patient)

We propose and evaluate two related mechanistic models based on a set of coupled ordinary differential equations (ODEs), either (a) deterministic or (b) stochastic. The stochastic differential equation (SDE) model comprises a Langevin equation, which can be converted into a Fokker-Planck partial differential equation (PDE) for the evolution of the probability distribution for pain gardiner1985handbook (). This allows for prediction of both the expected pain level for a patient at any point in the future and an assessment of the confidence in (and a confidence interval for) that prediction.

Mathematically, the deterministic mechanistic model we propose is the following, for a single patient:


where is the patient pain level (on a scale of 1–10), is the pain relaxation rate without drugs, is the marginal effect on the pain relaxation rate due to drug (), is the unmitigated pain level (i.e. without drug intervention), is the amount of standard drug doses within the patient, is the elimination rate of drug within the patient, are the drug dosage times, and is the number of doses of drug taken. represents the Dirac delta function. Note that the parameters and variables will in general need to be indexed with distinct values for each patient in a population, though we omit those indices here for clarity. For convenience and clarity we also here omit “hatted” notation (e.g., ) sometimes used for model predictions.

Tables 2 and 3 summarize the meanings of model variables and parameters, respectively.

Variable Meaning Units
Instantaneous pain level on 0–10 scale [pain]

Concentration of drug 1 (long-acting opioid) in the body [standard doses]

Concentration of drug 2 (short-acting opioid) in the body [standard doses]

Concentration of drug 3 (non-opioid) in the body [standard doses]

Instantaneous probability distribution of pain level [probability]
Table 2: Variables in mechanistic models.
Parameter Meaning Units
unmitigated pain level [pain]

rate of decrease of pain in the absence of drugs or acute sources of pain []

effect of drug 1 (long-acting opioid) on pain relaxation rate []

effect of drug 2 (short-acting opiod) on pain relaxation rate []

effect of drug 3 (non-opioid) on pain relaxation rate []

rate of decay of drug 1 (long-acting opioid) in body due to metabolism []

rate of decay of drug 2 (short-acting opioid) in body due to metabolism []

rate of decay of drug 3 (non-opioid) in body due to metabolism []

amplitude of intrinsic variability in human subjective pain reports

number of standard drug doses taken [count]

drug dose times (indexed by )
Table 3: Parameters in mechanistic models.

In this very simple model for pain dynamics (1), pain is expected to relax at rate to unmitigated level set by aggravating factors (like sickle cell disease) in the absence of intervention through opioids (drugs 1 and 2) or non-opioids (drug 3). When drugs are present in the patient’s body, pain drops at a faster rate and the short-term equilibrium pain level (not the unmitigated pain level ) is reduced. Note that we treat all parameters as constant over the time period of interest, which we take to be two weeks (based on clinical heuristic experience).

In the model for drug concentrations, medication in the body is assumed to be metabolized at a constant rate. Rates can be determined from existing substantiated pharmacokinetic models (e.g., yang2006 (); poulin2002 ()); Dirac delta function onset of medication serum concentration is a good approximation to the fast rise typical of the medications under consideration. See Fig. 2 for a sample deterministic model output.

Note that we deliberately chose to employ an extremely simple conceptual model for pain dynamics. More sophisticated versions might be developed to incorporate higher order dynamics for , or to include nonlinear or nonautonomous effects (e.g., allowing for explicit parameter variation with time of day or year), but currently available data are insufficient to constrain a model of greater complexity.

The stochastic differential (Langevin) equation version of our mechanistic model is as follows:


where a hypothesis of uncorrelated additive white noise has been made. From this we derive the Fokker-Planck equation for the probability distribution of pain over time :


Absent any pain medication, this Fokker-Planck equation implies the steady-state pain distribution


a Gaussian distribution with mean and standard deviation . See Figure 4 for a sample stochastic model output.

Figure 4: Sample output from stochastic differential equation model (2). Red thick line: theoretical mean pain; red thin lines: one theoretical standard deviation; black thick line: mean of pain distribution in ensemble of 100 stochastic simulations; blue thin lines: one standard deviation in ensemble of 100 stochastic simulations; blue dashed line: drug 1 dose in bloodstream. Spikes occur when patient takes recommended dosage.

ii.2.2 Statistical component (for all patients)

In order to account for the variation among patients and improve prediction of the unmitigated pain level, we associate patient characteristics and history with the unmitigated pain level (an -dimensional vector with corresponding to the -th patient’s unmitigated pain level) using a linear model.

Let be an design matrix containing the covariates of patients (i.e., patient characteristics). We write , with corresponding to the -dimensional covariates of patient . Then we formulate the relationship between between and the predictors as:


where is a -dimensional coefficient vector, and is an -dimensional vector of zero-mean random errors. When is small, the estimate for is obtained using the ordinary least squares procedure: , where denotes the norm. Then the unmitigated pain level is updated by , .

Since the unmitigated pain levels are not observable from patient pain reports, the initial ’s are independently sampled from a uniform distribution between and , i.e. . After using as initial values to fit the mechanistic model, the resulting estimated will be updated by the linear model (5) as , which will then be used as initial values in the next round of fitting of the mechanistic model. See Section II.3 for more detail on the hybridization of the statistical component with the mechanistic component.

Given a high-dimensional set of patient characteristics, we need to select a subset of patient characteristics that are significantly associated with by minimizing the penalized loss function. In this study, we select patient characteristics using the LASSO (Least Absolute Shrinkage and Selection Operator) tibshirani1996regression (), by minimizing the penalized loss function with respect to . The penalty parameter is determined using 5-fold cross-validation. The selected features are then used to fit the linear model (5) by ordinary least squares.

If time-varying unmitigated pain levels and time-varying covariates are present, the regression model (5) can be extended to the linear mixed model diggle2002analysis (); fitzmaurice2012applied (): where is an design matrix for random effect factors and is a vector of random effects. Patient characteristics can be selected by maximizing the penalized log-likelihood: groll2014variable (). Such an extension of the proposed hybrid model to allow for time-varying unmitigated pain levels and covariates will be considered in a future study with more data available.

ii.3 Model fitting

We fit our model to real patient data by minimizing the residual sum of squares between model predictions and patient reports provided within the first two weeks of reporting. We expect that the assumption of constant model parameters breaks down after approximately two weeks (clinical heuristics). Minimization over parameters , , and was done via the Nelder-Mead simplex algorithm nelder1965simplex (). Parameter was fixed at corresponding to a pain equilibration half-life time scale of 30 minutes in the absence of medication. If a patient did not take all three classes of drugs, the model and fitting only included the consumed drugs.

We initialize the parameter optimization in mechanistic models (one per patient) with random values during a first iteration, then we feed the optimization output into the statistical model (for all patients). Once the statistical model is run, it results in a new set of parameter estimates that can then be employed as initial parameter seeds for a second round of optimization in mechanistic models (to minimize the residual sum of squares). Proceeding iteratively in this fashion (see Fig. 3), we find convergence to a consistent set of parameters for each patient (details below).

ii.4 Method verification

Before applying our hybrid model to real-world patient data, we verify the soundness of the approach with synthetic data constructed to resemble real-world data, but generated by the model itself with high sampling frequency. The synthetic data used for verification of the method are generated directly from the mechanistic model with an assumed parameter set generated in the following way: unmitigated pain , initial pain level , and drug parameters . Each patient reports pain every 1/2 hour for 336 hours (two weeks). At each report time, the probability of the patient taking a particular drug (among three drug classes) is 1/16; in other words, the patients took each drug on average every 8 hours. White noise of magnitude 1 is added to each pain report.

As an illustration using real patient drug times (specifically those of Patient A3), we create synthetic data generated using and : see Figure 5. When the initial parameter search is seeded with random parameter values, the mechanistic model fit can lead to convergence to either the “true” optimum or to other “spurious” optima with incorrect values of , , and .

In this illustrative example, the method converges to and . The relevant root-mean-square (RMS) error is 1.01; this is close to the lowest possible expected error given the unit magnitude noise added to the synthetic data. This numerical experiment shows that the mechanistic model fitting method can converge even in the presence of significant amounts of noise. However, with only the mechanistic model it can be quite difficult to find a good set of initial parameter seeds111The seeding problem becomes exponentially harder as the dimension of the parameter space increases.: that is one motivation for introducing the statistical model.

Figure 5: Model fitting demonstration for densely reported noisy synthetic data. Upper panel: hypothetical densely-reported patient pain (black circles) and model fit (red solid line); red shading indicates model fit plus/minus one standard deviation. Lower panel: long-acting opioid (red solid line) and short-acting opioid (blue dashed line) medication concentration in patient bloodstream.

To test our hybrid method using both the mechanistic model for fitting and the statistical model for parameter estimation, we create a synthetic patient database of 39 patients as described above. We then iterate rounds of fitting between mechanistic and statistical models, starting with uniform random guesses for all patient parameters . Figure 6 demonstrates how the parameter converges to a value with small error after just a few iterations steps, even in the presence of significant noise. In order to evaluate the performance of the model on new data, we use the hold-out validation method by splitting the dataset into a training set (first week) and a test set (second week). Model fit error and hold-out validation error, as well as other parameters values, converge similarly: see Figure 7.

Figure 6: Hybrid model fitting demonstration for ensemble of densely reported noisy synthetic data. For an ensemble of 39 synthetic patient data sets, the average absolute error in gradually decreases. Iteration 0 indicates one fit to the mechanistic model alone. Subsequent iterations indicate the number of hybrid model (statistical + mechanistic) fits.
Figure 7: Hybrid model fitting demonstration for ensemble of densely reported noisy synthetic data. For an ensemble of 39 synthetic patient data sets, the average root-mean-squared (RMS) error in patient pain levels gradually decreases. Iteration 0 indicates one fit to the mechanistic model alone. Subsequent iterations indicate the number of hybrid model (statistical+mechanistic) fits. Training error (or fit error) is on the left; test error (or validation error) is on the right. Due to the additive white noise of magnitude 1, the smallest testing or training error we could expect is 1.

Iii Results

iii.1 General results

One key result is that our model for chronic pain does indeed have some predictive value (see Fig. 8). This is an improvement over the current state of the art, since no other predictive model exists of which we are aware. Furthermore, fitted parameter values correlate significantly with patient characteristics, suggesting that meaningful information is captured by this minimal plausible model. It may be possible to motivate new clinical insight on the basis of the observed correlations, perhaps leading to differential treatment of SCD sufferers with differing characteristics.

Figure 8: Hybrid model fitting on real patient data. For an ensemble of 39 real patient data sets, the average root-mean-squared (RMS) error in patient pain levels gradually decreases. Iteration 0 indicates one fit to the mechanistic model alone. Subsequent iterations indicate the number of hybrid model (statistical+mechanistic) fits. Training error (or fit error) is on the left; test error (or validation error) is on the right.

iii.2 Statistical results

We use the following baseline patient characteristics to predict the unmitigated pain levels in the statistical modeling step: age, gender, SCD disease type, hydroxyurea use, folic acid vitamin use, long-acting opioid use, short-acting opioid use, and non-opioid use. We explore the marginal effects of these characteristics and their possible pairwise two-way interactions using the LASSO. The model (5) can be extended to include time-varying covariates such as temperature, weather, patient’s walking/social activities, and patient’s mood at time , once these data become available in a future study.

The statistical model that resulted from the LASSO variable selection is given by


where , , and is the indicator function.

Table 4 summarizes the results from one round of fitting of the regression model (6). Adjusting for the effect of other terms in the regression model, SCD disease type of SBThal or So-Ara (with coefficient ), non-opioid use (with coefficient ), and the interaction term between SCD disease type of SBThal or So-Ara and age (with coefficient ) are important predictors of the unmitigated pain levels at the significant level of 0.05. Using non-opioid medication is associated with decreased unmitigated pain levels. Unmitigated pain levels increase with patients’ age for SBThal or So-Ara patients.

Variable Estimate Std Err T-value P-value
Intercept 7.646 1.228 6.228 0.000 ***
HgbSC -1.566 0.890 -1.761 0.088
SBThal or So-Ara -5.479 2.332 -2.349 0.025 *
Age at baseline 18 0.001 0.034 0.290 0.773
Hydroxyurea user -1.205 0.839 -1.437 0.160
Non-opioid user -2.523 0.842 -2.995 0.005 **
(SBThal or So-Ara) (Age at baseline -18) 0.241 0.010 2.419 0.021 *
Table 4: Result of the prediction model of the unmitigated pain using the linear regression model. Significance

iii.3 Mechanistic model validation

With such sparse data and up to four fitting parameters, one may worry that the model (1) is being overfitted. To test this concern, we propose 6 related alternative models with fewer fitting parameters, and we compare cross-validation error and Akaike information criterion (AIC) among the models. See Table 5 for model descriptions. Neither measure selected a best-fit model across all patients, but none of these simple models is overfitting the data. See Figure 9 for AIC results and Figure 10 for cross validation results.

Model name Description Fitting parameters
Full model Include all drugs taken (model (1)) up to 4

No drugs
Include no drug dosing information 1

Merge drugs
Combine all drugs into one drug class with same response up to 2

LA only
Include only long-acting opioid doses up to 2

SA only
Include only short-acting opioid doses up to 2

NO only
Include only non-opioid doses up to 2

Include drug class only if drug is taken at least times* up to 4
Table 5: Mechanistic model variations. Fitting parameters include unmitigated pain level and drug response parameters for all drugs consumed. Therefore some patients have fewer fitting parameters than listed if they consumed fewer than three types of drugs. *In our tests, was the drug dose threshold.
Figure 9: Akaike information criterion (AIC) for alternative models listed in Table 5. Most models perform equally well; among patients with differing model performance, there exists no clear ‘best’ model for all patients.
Figure 10: Two-fold cross validation testing error for alternative models listed in Table 5. For every patient, each model was independently fitted to the first and second half of the time series pain report data (training). Then the fitted models were used to test the other half of the data. This figure shows the average root-mean-square testing error for the two tests. Most models perform equally well; among patients with differing model performance, there exists no clear ‘best’ model for all patients. Note that patient #36 did not have enough data to fit any models, so zero error is misleading.

Iv Biased pain reporting

Perhaps the most significant limitation of our model lies in a potential bias in our data set. Patients typically report pain levels when taking medication, but many of them only take medication when pain levels rise. Thus we suspect a selection bias of unknown significance, causing higher pain levels to be reported at a disproportionately high rate. To test this concern, we compare the unbiased model (1) with a similar model incorporating biased pain reporting.

Suppose the probability density function of pain at a particular time is a normal distribution with mean and variance :


Integrating model (1) gives the expected pain value at any point in time.

If higher pain is disproportionately reported through the mobile health application, then we will be much more likely to see higher pain levels from this normal distribution. As a first approximation, we assume the reporting bias is linear:


where normalizes the distribution, tune the probabilities of reporting a pain value , and the Heaviside function prevents negative pain values. Figure 11 shows both real and reported pain distributions at a particular time.

Figure 11: Probability density distributions for unbiased (solid) and biased (dashed) pain reporting at a particular time . The standard deviation (here, ) has been exaggerated for illustrative purposes. In real data, the typical standard deviation is around .

We need a way to connect these distributions because we want to control real pain described by (7), but the patient only provides data from the reported distribution (8). In other words, real pain is important but invisible, and reported pain is unimportant but visible. One way to connect the distributions is through their means and variances222Note that the means and variances also change in time. We omit time dependence for notational clarity.:


Assuming is approximately constant in time333It is not possible to verify this with data because patients only report one pain value at a particular time. However, a Kolmogorov-Smirnov normality test on residuals over the first two weeks of data rejects normality () for only 2 of the 39 patients., we can also estimate using the definition of variance:


where is the number of pain reports, is the th reported pain value, and is the expected reported pain given distribution (8) at time .

Given a proposed model for real pain, we can solve this system of three equations for the three unknowns: , , and . We can then compute the likelihood of the reported pain using


Because we can compute the likelihood of the supplied data given any proposed model for real pain, we can tune the model parameters to maximize likelihood (technically, we minimize the negative logarithm of likelihood). This results in a best-fit model under the assumption of biased pain reporting. We compare the best-fits of the model under both biased and unbiased reporting assumptions, and find that neither model is a better fit for most patients. See Figure 12.

Figure 12: Akaike information criterion for unbiased and biased pain reporting. Because both models have an equal number of fitting parameters, AIC is a proxy for model likelihood (lower AIC implies higher likelihood). Again, it is not clear that one model performs universally better than the other. Note that missing biased reporting model fits indicate that the fitting algorithm did not converge ().

V Pain and medication optimization

A key goal of the modeling of human pain dynamics is to develop predictions that allow optimized treatment: both pain and medication use should be minimized. Excess medication carries particular long-term risks for chronic pain sufferers bannwarth1999risk (); gatchel2001biopsychosocial (); savage1999opioid (); brookoff2000chronic (), but pain mitigation is also a primary goal of SCD treatment. How can these contradictory objectives be balanced?

Our model allows us to forecast the probability distribution of pain for a patient at a point in the near future, given past data and future drug dosage protocol. This information may be useful to a physician, allowing him or her to make an optimized, data-driven decision balancing medication and pain for the patient in real time.

We propose several tools that may be useful to a physician. First, we find the optimal drug timing given that a certain amount of each drug will be taken over a certain time period (say, within 24 hours). For instance, if the patient will take two short-acting opioids and one long-acting opioid within 24 hours, then the algorithm will offer the best times to take those three drug doses in order to minimize the expected average pain. We provide the physician with the expected optimal average pain for all drug combinations up to a certain maximum number of safe drug doses. See Figure 13.

Figure 13: Example expected pain given optimal drug dosage protocol (in this case, for Patient A3). For each set of drug dosage protocols, from no drugs to 4 doses of each drug, a physician can see the expected average pain over a certain time period. Given a patient’s maximum acceptable pain level, the physician can select the best compromise between drug doses and expected pain. In this case, the physician may tell the patient to take no long-acting (LA) drugs but take 3-4 short-acting (SA) drugs. Alternatively, the physician might tell the patient to take one LA and 1-2 doses of SA medication. The timing of those LA and SA drugs is provided by the optimization algorithm.

Second, we select the best drug dosing protocol (both number of drugs and dose timings) given an objective function balancing pain and medication. There are an unlimited number of possible objective functions that balance pain and medication, but we propose the following:


where is the average expected pain; are the number of drug doses of each type; are the weights of pain and drugs; are the maximum safe levels of pain and drugs; and tune the steepness of the objective function near those dangerous levels.

See Figure 14 for the contributions of the pain and one drug component to the objective function . The contribution to the objective function is zero if no pain exists or if no drugs are taken. As pain or drug doses approach dangerous levels, the contribution to the objective function blows up. After a physician has made sufficient recommendations to a patient, a machine learning algorithm could select the weight parameters for each physician/patient pair. At that point, the algorithm could propose the optimal dosing protocol without much effort on the physician’s part. See Figure 15.

Figure 14: Contribution of the pain and one drug component to the objective function (13). The contribution to the objective function is 0 if no pain exists or if no drugs are taken. As pain or drug doses approach dangerous levels, the contribution to the objective function blows up.
Figure 15: Example patient intervention recommendation (in this case, for Patient A3). For a set of personalized optimization parameters (selected by a physician or machine learning algorithm), the optimal drug timing minimizes the objective function (13) for each number of drug doses per 24-hour period. In this case, the patient is advised to take two standard doses of the long-acting (LA) opioid and two standard doses of the short-acting (SA) opioid, indicated by the red box.

Vi Discussion and Limitations

vi.1 Reflection on hybrid modeling

Statistical models and mechanistic models have both been successfully applied to various aspects of human behavior. The inference of “black box” statistical models from empirical data has the advantage that it obviates the need for a-priori knowledge of system dynamics. However, mechanistic models (sometimes referred to as “white box” or “clear box”) can easily incorporate such knowledge when available.

Perhaps because of the often distinct educational backgrounds of practitioners or distinct typical applications, statistical and mechanistic approaches are not frequently combined in addressing a single problem. Compared with our work, the most similar hybrid modeling idea was developed by Sheiner and colleagues in the field of pharmacokinetics, where they proposed models to estimate population characteristics of pharmacokinetic parameters sheiner1977estimation (); sheiner1980evaluation (); mandema1992building (). In their work, the pharmacokinetic models (i.e., mechanistic models) are well established, and the novelty and focus was the introduction of statistical models for pharmacokinetic parameter estimation. On the contrary, in our study the mechanistic model is not known before but developed by us based on clinical knowledge and reasonable assumptions, and our focus is the prediction of pain levels rather than parameter estimation.

Other attempts based on the hybrid modeling idea in the scientific literature have appeared in the context of neural networks (e.g., psichogios1992hybrid (); thompson1994modeling (); su2014integrating ()) and chemical engineering (e.g., thompson1994modeling (); schubert1994hybrid (); duarte2003hybrid ()), where they largely played a computational rather than analytical role. Some attempts have also been made with medical applications: Rosenberg et al. (rosenberg2007using ()) and Adams et al. (adams2007estimation ()) developed a model by combining a dynamical systems approach with a statistical model to predict a patient’s CD4 cell counts and HIV viral load over time in an HIV study. Timms et al. (timms2014dynamical ()) proposed a dynamical systems approach using ODEs to improve self-regulation in a smoking cessation study. Reinforcement learning techniques such as Q-learning (e.g., jaimes2014stress ()) also share some commonalities with the hybrid approach.

In this work we make our own attempt at a novel incorporation of statistical inference together with mechanistic dynamical systems modeling to produce a hybrid mathematical model for predicting and explaining human behavior. We apply the new approach specifically to the problem of predicting the dynamics of subjective pain in a population of individuals suffering from sickle cell disease. The rationale behind our method development is that we have prior knowledge of pain trajectories with medication, making the problem suitable for mechanistic modeling; meanwhile, we do not know the relationship between patient health characteristics and pain levels, so we would like to investigate this using a statistical model.

vi.2 Limitations

The hybrid dynamical systems/statistical approach appears to have great potential. The low frequency of pain reporting currently limits its usefulness, but future addition of high-frequency pain correlates like blood pressures, heart rate, activity level, etc., via wearable medical devices (e.g. the “Fitbit”) may drastically improve on that. Furthermore, application of similar methods to more data-rich forecasting problems (e.g. insulin levels) may also expand the utility of our work.

Another important limitation to our current model lies in the mechanistic component. We presented here what we considered to be the simplest plausible model: pain fluctuates about an “unmitigated” equilibrium , and medication reduces pain below that level, but pain returns as medication is metabolized and removed from the bloodstream. This simple model cannot capture long-term changes in the unmitigated pain level, and hence its forecast validity is likely limited to short time scales (days to weeks).

Society is clearly moving in the direction of an overwhelming amount of medical data. It is imperative that we learn to take advantage of this information to improve patient treatments beyond the traditional standard of care. The approach we report here not only addresses the specific challenge of chronic pain mitigation in SCD patients, but also provides a testbed for new ways of dealing with big, ever-growing data sets in real time.

Vii Conclusions

We have successfully demonstrated the hybrid application of statistical and mechanistic mathematical modeling with application to understanding the dynamics of subjective human pain. Our model explains real-world data on human pain and can generate predictions of future pain dynamics.

We expect that similar methods could be used to incorporate disease-specific knowledge and modeling with statistical inference in a variety of medical applications. Given the coming deluge of data from wearables (including clinical trial NCT02895841 already underway) and mobile health applications, there is a clear need for new mathematical methods to take advantage of the opportunity for personalizable, data-driven medical treatments.

Viii Data Accessibility

Data and code available at Northwestern University ARCH repository: https://doi.org/10.21985/N24S3T NUDataSoftware (). Real patient data is not included to protect patient privacy, but all methods may be tested on synthetic data.

Ix Funding

The authors gratefully acknowledge NSF support through grants No. DMS-1557727 and No. DGE-1324585.


  • (1) S. H. Strogatz, Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. Westview press, 2014.
  • (2) D. A. Freedman, Statistical models: theory and practice. cambridge university press, 2009.
  • (3) O. S. Platt, B. D. Thorington, D. J. Brambilla, P. F. Milner, W. F. Rosse, E. Vichinsky, and T. R. Kinney, “Pain in sickle cell disease: rates and risk factors,” New England Journal of Medicine, vol. 325, no. 1, pp. 11–16, 1991.
  • (4) Y. Granovsky, M. Granot, R.-R. Nir, and D. Yarnitsky, “Objective correlate of subjective pain perception by contact heat-evoked potentials,” The Journal of Pain, vol. 9, no. 1, pp. 53–63, 2008.
  • (5) A. Hughes, A. Macleod, J. Growcott, and I. Thomas, “Assessment of the reproducibility of intradermal administration of capsaicin as a model for inducing human pain,” Pain, vol. 99, no. 1, pp. 323–331, 2002.
  • (6) S. Stevens, A. Carton, and G. Shickman, “A scale of apparent intensity of electric shock.,” Journal of Experimental Psychology, vol. 56, no. 4, p. 328, 1958.
  • (7) C. R. Jonassaint, N. Shah, J. Jonassaint, and L. De Castro, “Usability and feasibility of an mhealth intervention for monitoring and managing pain symptoms in sickle cell disease: The sickle cell disease mobile application to record symptoms via technology (smart),” Hemoglobin, no. ahead-of-print, pp. 1–7, 2015.
  • (8) N. Shah, J. Jonassaint, and L. De Castro, “Patients welcome the sickle cell disease mobile application to record symptoms via technology (smart),” Hemoglobin, vol. 38, no. 2, pp. 99–103, 2014.
  • (9) D. T. McRuer and E. S. Krendel, “Mathematical models of human pilot behavior,” tech. rep., DTIC Document, 1974.
  • (10) U. G. Fors, L. G. Edwall, and G. A. Haegerstam, “The ability of a mathematical model to evaluate the effects of two pain modulating procedures on pulpal pain in man,” Pain, vol. 33, no. 2, pp. 253–264, 1988.
  • (11) N. Britton, S. Skevington, and M. Chaplain, “Mathematical modelling of acute pain,” Journal of Biological Systems, vol. 3, no. 04, pp. 1119–1124, 1995.
  • (12) G. Stepan, “Delay effects in the human sensory system during balancing,” Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, vol. 367, no. 1891, pp. 1195–1212, 2009.
  • (13) C. W. Gardiner et al., Handbook of stochastic methods, vol. 3. Springer Berlin, 1985.
  • (14) F. Yang, X. Tong, D. G. McCarver, R. N. Hines, and D. A. Beard, “Population-based analysis of methadone distribution and metabolism using an age-dependent physiologically based pharmacokinetic model,” Journal of pharmacokinetics and pharmacodynamics, vol. 33, no. 4, pp. 485–518, 2006.
  • (15) P. Poulin and F.-P. Theil, “Prediction of pharmacokinetics prior to in vivo studies. ii. generic physiologically based pharmacokinetic models of drug disposition,” Journal of pharmaceutical sciences, vol. 91, no. 5, pp. 1358–1370, 2002.
  • (16) R. Tibshirani, “Regression shrinkage and selection via the lasso,” Journal of the Royal Statistical Society. Series B (Methodological), pp. 267–288, 1996.
  • (17) P. Diggle, Analysis of longitudinal data. Oxford University Press, 2002.
  • (18) G. M. Fitzmaurice, N. M. Laird, and J. H. Ware, Applied longitudinal analysis, vol. 998. John Wiley & Sons, 2012.
  • (19) A. Groll and G. Tutz, “Variable selection for generalized linear mixed models by l 1-penalized estimation,” Statistics and Computing, vol. 24, no. 2, pp. 137–154, 2014.
  • (20) J. A. Nelder and R. Mead, “A simplex method for function minimization,” The computer journal, vol. 7, no. 4, pp. 308–313, 1965.
  • (21) B. Bannwarth, “Risk-benefit assessment of opioids in chronic noncancer pain,” Drug safety, vol. 21, no. 4, pp. 283–296, 1999.
  • (22) R. J. Gatchel, “A biopsychosocial overview of pretreatment screening of patients with pain,” The Clinical journal of pain, vol. 17, no. 3, pp. 192–199, 2001.
  • (23) S. Savage, “Opioid therapy of chronic pain: assessment of consequences,” Acta Anaesthesiologica Scandinavica, vol. 43, no. 9, pp. 909–917, 1999.
  • (24) D. Brookoff, “Chronic pain: 2. the case for opioids,” Hospital Practice, vol. 35, no. 9, pp. 69–84, 2000.
  • (25) L. B. Sheiner, B. Rosenberg, and V. V. Marathe, “Estimation of population characteristics of pharmacokinetic parameters from routine clinical data,” Journal of pharmacokinetics and biopharmaceutics, vol. 5, no. 5, pp. 445–479, 1977.
  • (26) L. B. Sheiner and S. L. Beal, “Evaluation of methods for estimating population pharmacokinetic parameters. i. michaelis-menten model: routine clinical pharmacokinetic data,” Journal of pharmacokinetics and biopharmaceutics, vol. 8, no. 6, pp. 553–571, 1980.
  • (27) J. W. Mandema, D. Verotta, and L. B. Sheiner, “Building population pharmacokineticpharmacodynamic models. i. models for covariate effects,” Journal of pharmacokinetics and biopharmaceutics, vol. 20, no. 5, pp. 511–528, 1992.
  • (28) D. C. Psichogios and L. H. Ungar, “A hybrid neural network-first principles approach to process modeling,” AIChE Journal, vol. 38, no. 10, pp. 1499–1511, 1992.
  • (29) M. L. Thompson and M. A. Kramer, “Modeling chemical processes using prior knowledge and neural networks,” AIChE Journal, vol. 40, no. 8, pp. 1328–1340, 1994.
  • (30) H.-T. Su, N. Bhat, P. Minderman, and T. McAvoy, “Integrating neural networks with first principles models for dynamic modeling,” in Dynamics and Control of Chemical Reactors, Distillation Columns and Batch Processes (DYCORD+’92): Selected Papers from the 3rd IFAC Symposium, Maryland, USA, 26-29 April 1992, p. 327, Elsevier, 2014.
  • (31) J. Schubert, R. Simutis, M. Dors, I. Havlík, and A. Lübbert, “Hybrid modelling of yeast production processes–combination of a priori knowledge on different levels of sophistication,” Chemical engineering & technology, vol. 17, no. 1, pp. 10–20, 1994.
  • (32) B. P. Duarte and P. M. Saraiva, “Hybrid models combining mechanistic models with adaptive regression splines and local stepwise regression,” Industrial & engineering chemistry research, vol. 42, no. 1, pp. 99–107, 2003.
  • (33) E. S. Rosenberg, M. Davidian, and H. T. Banks, “Using mathematical modeling and control to develop structured treatment interruption strategies for hiv infection,” Drug and alcohol dependence, vol. 88, pp. S41–S51, 2007.
  • (34) B. Adams, H. Banks, M. Davidian, and E. Rosenberg, “Estimation and prediction with hiv-treatment interruption data,” Bulletin of mathematical biology, vol. 69, no. 2, pp. 563–584, 2007.
  • (35) K. P. Timms, D. E. Rivera, L. M. Collins, and M. E. Piper, “A dynamical systems approach to understanding self-regulation in smoking cessation behavior change,” nicotine & tobacco research, vol. 16, no. Suppl 2, pp. S159–S168, 2014.
  • (36) L. G. Jaimes, M. Llofriu, and A. Raij, “A stress-free life: just-in-time interventions for stress via real-time forecasting and intervention adaptation,” in Proceedings of the 9th International Conference on Body Area Networks, pp. 197–203, ICST (Institute for Computer Sciences, Social-Informatics and Telecommunications Engineering), 2014.
  • (37) S. M. Clifton, C. Kang, J. J. Li, Q. Long, N. Shah, and D. M. Abrams, “Data from: Hybrid statistical and mechanistic mathematical model guides mobile health intervention for chronic pain.”
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description