As one might see in Fig. \ref{198684}, after we constrained most of the parameters to reasonable values, it became more difficult to achieve a satisfactory fit. In particular, while it wasn't difficult to generate plausible scenarios, it proved difficult to replicate the explosive growth patterns in the death rate or the prescribed population.
One possible reason for this might be the difficulty in optimizing a complex model: the term for the addicted individuals in the model, \(A\), has a complex relationship with every other group, some nonlinearly, all of which constrain the system. This means it's difficult to get all the factors to align in order to get the desired relationship. While using a multi-start method improved our results, it's possible that we are still not achieving a global minimum. One reason for the model fit we see in the case of the prescribed population may be inaccuracies in our data or response variable: we currently predict the prescribed population over time by combining data for the number of prescriptions from one data source with the estimated number of prescriptions per person from another data source. Data that directly reports the number of prescribed users over time would be helpful in achieving better model results. Finally, the fact that we are trying to optimize the fit to multiple time series at the same time presents an additional constraint to this method.

Conclusions

In this paper, we discuss a mathematical approach for analyzing the opioid epidemic, through extending and improving the model from Battista et al. (2017).  We sought first to clarify the terms and assumptions in this model. Through applying an additional assumption to the nonlinear terms, and removing an assumption for which we found no evidence, we were able to reduce the total number of parameters in the model and leave only terms that are easily interpretable.
We present a method for fitting our model to real world data sources, which has not been attempted in cited work. We encountered many challenges in finding parameters that simulate the explosive growth of opioid prescriptions and deaths that we see in the data. While continued iterations of the model and its assumptions may improve this fit, we also believe there is significant possibility for inaccuracy in our data which may be hindering the result. Future research may involve further investigation into assessing which aspects of the model or data lead to the fit we obtained. While we employed multi-start optimization methods in hopes of finding a global minimum, we suggest applying other methods for global optimization in the future as well. 
Through global sensitivity analysis, we found that the parameters most sensitive to the endemic addicted and recovery populations were  \(\zeta\), the rate at which addicts enter treatment, and \(\delta\), the rate at which recovering addicts successfully complete treatment. We found that model parameters related to new addictions, such as \(\beta_1\)\(\gamma\) and \(\alpha\), are much less impactful. This suggests that encouraging addicted individuals to seek treatment and working to increase its effectiveness are of utmost importance in reducing the total number of people affected by this crisis. While preventing new addictions is of course helpful, our result underscores the fact that we currently a large population of individuals stuck in addiction. Interventions, then, must directly address rehabilitating that population - there is no other pathway through which these individuals can move back to the susceptible group and be free of addiction.