Hybrid statistical and mechanistic mathematical model guides mobile health intervention for chronic pain
Abstract
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 realtime 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 realtime adaptive treatment recommendations for pain. Furthermore, current stateoftheart methods have difficulty in handling continuoustime 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 apriori 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 patientreported pain levels vary over time based on medication dosage information and other patient characteristics.
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 30day hospital reutilization rates are alarmingly high platt1991pain (). While factors influencing these high reutilization rates are poorly understood, close followup and continued use of pain medications has been shown to decrease rehospitalization rates. Mobile technology has become an integral part of health care management, and our recently selfdeveloped 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 nonopioid 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 010 scale) in real time jonassaint2015usability (); shah2014patients ().
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  
N  
Institution  
A  14  
B  17  
C  8  
Gender  
Male  16  
Female  23  
Age at baseline (years)  
1834  24  
15  
SCD disease type  
Hemoglobin SC  8  
Hemoglobin SS  22  
Hemoglobin SB (Beta) Thalassemia  5  
BetaZero Thalassemia  3  
SOAra  1  
Hydroxyurea user  27  
Folic acid vitamin user  26  
Longacting opioid user  29  
Shortacting opioid user  35  
Nonopioid 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 )  
Withinpatient 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 longacting opioid doses (first 2 weeks)  6.0  ( 0.0, 35.0 )  
Number of shortacting opioid doses (first 2 weeks)  7.2  ( 0.0, 35.0 )  
Number of nonopioid doses (first 2 weeks)  2.1  ( 0.0, 12.0 )  

ii.2 Predictive model
In order to develop a hybrid model that incorporates both a mechanistic apriori knowledgedriven component and a statistical datadriven 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 patientspecific 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.
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 FokkerPlanck 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:
(1)  
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.
Variable  Meaning  Units 

Instantaneous pain level on 0–10 scale  [pain]  

Concentration of drug 1 (longacting opioid) in the body  [standard doses] 

Concentration of drug 2 (shortacting opioid) in the body  [standard doses] 

Concentration of drug 3 (nonopioid) in the body  [standard doses] 

Instantaneous probability distribution of pain level  [probability] 
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 (longacting opioid) on pain relaxation rate  [] 

effect of drug 2 (shortacting opiod) on pain relaxation rate  [] 

effect of drug 3 (nonopioid) on pain relaxation rate  [] 

rate of decay of drug 1 (longacting opioid) in body due to metabolism  [] 

rate of decay of drug 2 (shortacting opioid) in body due to metabolism  [] 

rate of decay of drug 3 (nonopioid) 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 ) 
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 nonopioids (drug 3). When drugs are present in the patient’s body, pain drops at a faster rate and the shortterm 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:
(2)  
where a hypothesis of uncorrelated additive white noise has been made. From this we derive the FokkerPlanck equation for the probability distribution of pain over time :
(3) 
Absent any pain medication, this FokkerPlanck equation implies the steadystate pain distribution
(4) 
a Gaussian distribution with mean and standard deviation . See Figure 4 for a sample stochastic model output.
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:
(5) 
where is a dimensional coefficient vector, and is an dimensional vector of zeromean 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 highdimensional 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 5fold crossvalidation. The selected features are then used to fit the linear model (5) by ordinary least squares.
If timevarying unmitigated pain levels and timevarying 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 loglikelihood: groll2014variable (). Such an extension of the proposed hybrid model to allow for timevarying 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 NelderMead simplex algorithm nelder1965simplex (). Parameter was fixed at corresponding to a pain equilibration halflife 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 realworld patient data, we verify the soundness of the approach with synthetic data constructed to resemble realworld 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 rootmeansquare (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 seeds^{1}^{1}1The seeding problem becomes exponentially harder as the dimension of the parameter space increases.: that is one motivation for introducing the statistical model.
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 holdout validation method by splitting the dataset into a training set (first week) and a test set (second week). Model fit error and holdout validation error, as well as other parameters values, converge similarly: see Figure 7.
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.
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, longacting opioid use, shortacting opioid use, and nonopioid use. We explore the marginal effects of these characteristics and their possible pairwise twoway interactions using the LASSO. The model (5) can be extended to include timevarying 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
(6) 
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 SoAra (with coefficient ), nonopioid use (with coefficient ), and the interaction term between SCD disease type of SBThal or SoAra and age (with coefficient ) are important predictors of the unmitigated pain levels at the significant level of 0.05. Using nonopioid medication is associated with decreased unmitigated pain levels. Unmitigated pain levels increase with patients’ age for SBThal or SoAra patients.
Variable  Estimate  Std Err  Tvalue  Pvalue  

Intercept  7.646  1.228  6.228  0.000  *** 
HgbSC  1.566  0.890  1.761  0.088  
SBThal or SoAra  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  
Nonopioid user  2.523  0.842  2.995  0.005  ** 
(SBThal or SoAra) (Age at baseline 18)  0.241  0.010  2.419  0.021  * 
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 crossvalidation error and Akaike information criterion (AIC) among the models. See Table 5 for model descriptions. Neither measure selected a bestfit 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 longacting opioid doses  up to 2 
SA only 
Include only shortacting opioid doses  up to 2 
NO only 
Include only nonopioid doses  up to 2 
Threshold 
Include drug class only if drug is taken at least times*  up to 4 
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 :
(7) 
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:
(8) 
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.
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 variances^{2}^{2}2Note that the means and variances also change in time. We omit time dependence for notational clarity.:
(9)  
(10) 
Assuming is approximately constant in time^{3}^{3}3It is not possible to verify this with data because patients only report one pain value at a particular time. However, a KolmogorovSmirnov 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:
(11) 
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
(12) 
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 bestfit model under the assumption of biased pain reporting. We compare the bestfits 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.
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 longterm 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, datadriven 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 shortacting opioids and one longacting 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.
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:
(13) 
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.
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 apriori 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 selfregulation in a smoking cessation study. Reinforcement learning techniques such as Qlearning (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 highfrequency 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 datarich 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 longterm 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, evergrowing 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 realworld data on human pain and can generate predictions of future pain dynamics.
We expect that similar methods could be used to incorporate diseasespecific 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, datadriven 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. DMS1557727 and No. DGE1324585.
References
 (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 heatevoked 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. aheadofprint, 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, “Populationbased analysis of methadone distribution and metabolism using an agedependent 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 1penalized 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, “Riskbenefit 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. michaelismenten 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 networkfirst 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, 2629 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 hivtreatment 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 selfregulation 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 stressfree life: justintime interventions for stress via realtime forecasting and intervention adaptation,” in Proceedings of the 9th International Conference on Body Area Networks, pp. 197–203, ICST (Institute for Computer Sciences, SocialInformatics 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.”