Methods for determining key components in a mathematical model for tumor–immune dynamics in multiple myeloma

Highlights

• We explore a mathematical model of multiple myeloma that was presented previously.

• We justify model parameter ranges and values using published data and calculations.

• We use global methods to determine parameters the model is most sensitive to.

• We determine which of those parameters could be estimated from data.

• We numerically explore the behavior of the model over time.

Abstract

In this work, we analyze a mathematical model we introduced previously for the dynamics of multiple myeloma and the immune system. We focus on four main aspects: (1) obtaining and justifying ranges and values for all parameters in the model; (2) determining a subset of parameters to which the model is most sensitive; (3) determining which parameters in this subset can be uniquely estimated given certain types of data; and (4) exploring the model numerically. Using global sensitivity analysis techniques, we found that the model is most sensitive to certain growth, loss, and efficacy parameters. This analysis provides the foundation for a future application of the model: prediction of optimal combination regimens in patients with multiple myeloma.

Introduction

Multiple myeloma (MM) is a cancer of plasma B cells. Although there are almost two dozen treatments approved in the US and others currently in clinical trials, many patients do not survive more than ten years (Kazandjian and Landgren, 2016). Regimens usually combine multiple therapies (two, three, or more drugs) simultaneously. Patients receive additional combination therapies if they do not respond to the initial combination, or when they relapse. Without direct comparison studies, it is difficult to know the best combinations and doses for treatment, and questions remain regarding both treatment choice and timing (Kazandjian and Landgren, 2016). A mathematical model that accurately captures important features of the disease and therapy dynamics can be used to test regimens in silico. It can also be used to calculate a regimen that is predicted to perform optimally.

A number of mathematical models for the progression of MM and its response to treatments have been developed previously. Most of these models have made use of the correlation between MM tumor burden and a protein that is shed by MM cells, M protein, that can be measured in the peripheral blood (Durie, Salmon, 1975, Salmon, Smith, 1970). Sullivan and Salmon (1972) published a simple mathematical model to characterize chemotherapy-induced tumor regression in patients with MM. Swan and Vincent applied optimal control in Swan and Vincent (1977) to predict an optimal chemotherapy dosing strategy for patients with MM. Hokanson et al. (1977) fit individual M-protein data with mathematical models of myeloma cell populations that were sensitive or resistant to chemotherapy. More recently, several groups have examined the disease progression in the presence of various treatments. Jonsson et al. (2015) published a model of patient M-protein levels in response to treatment with carfilzomib. Tang et al. (2016) published a model for patient M-protein levels in clinical trials of bortezomib-based chemotherapy. They predicted that rational combination treatments with decreased selection pressure on myeloma cells could lead to a longer remission period. Nanavati et al. (2017) published a semi-mechanistic protein production and signaling model with seven main compartments. The model was fit to aggregate in vivo xenograft data (digitized from the literature) for mice treated with vorinostat.

These prior studies have not directly examined the role of the immune system in MM disease dynamics, though there is a long history of tumor–immune dynamics models (cf. d’Onofrio, 2005, Kirschner, Panetta, 1998, Kuznetsov, Makalkin, Taylor, Perelson, 1994, Moore, Li, 2004, de Pillis, Radunskaya, Wiseman, 2005, Stepanova, 1980, de Vladar, González, 2004). One contribution of our work is the incorporation and careful parameterization of immune-disease interactions in a mathematical model for MM. This is important as several immunomodulatory drugs have been approved for use in patients with MM, and more are currently in clinical trials (Kazandjian and Landgren, 2016). A well-justified model of tumor–immune dynamics would support in silico exploration of regimens that include immunomodulatory therapies.

In Gallaher et al. (2018), we introduced a semi-mechanistic mathematical model of MM (tumor) and immune system dynamics in a hypothetical patient. We focused on the structure of the model, and how best to incorporate tumor–immune dynamics. The mathematical model consists of a dynamical system that tracks tumor and immune components in the peripheral blood of patients with MM. Although MM is a disease based in the bone marrow, myeloma cells overproduce a myeloma protein or M protein, which can circulate outside the bone marrow. Our model uses the level of M protein in the peripheral blood as a surrogate of tumor burden in patients with MM (Durie, Salmon, 1975, Salmon, Smith, 1970). The immune cells included in our model play important roles in disease control or progression, and can also be measured in the peripheral blood. In addition to justifying the structure of the full model, this work also considered a reduced version of the model that still captures key long-term dynamics, and analyzed equilibria and stability of this reduced model.

In this paper, we present methods to determine important components or interactions in the model. This information can be used to help decide which therapies to consider in combination regimens. We began with extensive literature searches to determine and justify parameter values and ranges to use in the model. For the ranges considered feasible for the parameters, we used sampling to perform global sensitivity analyses. We found eight parameters that have the largest effect on long-term levels of M protein. We refer to this subset as “sensitive parameters”. We numerically explored how small changes in the sensitive parameters can lead to qualitatively different types of model behavior. In particular, changes in values of the sensitive parameters can result in a switch between high and low tumor burden states.

With all except the sensitive parameters fixed, we used identifiability analysis to determine that the model is globally structurally identifiable. In other words, the eight sensitive parameters could be uniquely estimated if we had continuous data from all four populations in our model. Although we did not find time-series data in the literature, we did find steady-state patient values for the four tumor and immune populations included in our model (Greipp, Miguel, Durie, Crowley, Barlogie, Bladé, Boccadoro, Child, Avet-Loiseau, Kyle, Lahuerta, Ludwig, Morgan, Powles, Shimizu, Shustik, Sonneveld, Tosi, Turesson, Westin, 2005, Pessoa de Magalhães, Vidriales, Paiva, Fernandez-Gimenez, García-Sanz, Mateos, Gutierrez, Lecrevisse, Blanco, Hernández, de las Heras, Martinez-Lopez, Roig, Costa, Ocio, Perez-Andres, Maiolino, Nucci, Rubia, Lahuerta, San-Miguel, Orfao, 2013, Tang, Zhao, van de Velde, Tross, Mitsiades, Viselli, Neuwirth, Esseltine, Anderson, Ghobrial, Miguel, Richardson, Tomasson, Michor, 2016). Fitting to this small amount of data, we were able to explore relationships between the sensitive parameters.

We propose a model for tumor–immune dynamics of multiple myeloma that can be used in future work. The model we propose was developed in Gallaher et al. (2018), has parameter values/ranges we based on literature searches, and has fixed values for all but the eight most sensitive parameters found in this current work. The analysis and findings in this and previous work (Gallaher et al., 2018) add to our overall confidence in the model. If individual-level clinical data become available for each of the four populations in the model, the eight sensitive parameters could be estimated. In silico exploration and optimization could then identify regimens predicted to have best outcomes. Such regimens could be tested either preclinically or clinically and compared to others to evaluate outcomes and calibrate the results.