Model based analysis of the heterogeneity in the tumour size dynamics differentiates vemurafenib, dabrafenib and trametinib in metastatic melanoma

Purpose Explore the heterogeneity in dynamics of tumour response to vemurafenib, dabrafenib and trametinib using routinely collected clinical trial imaging data. Methods Time-series imaging data from the phase III studies of vemurafenib, dabrafenib and trametinib were collected through a data repository. A mathematical model based on basic mechanisms of tumour growth was placed within a statistical modelling framework to analyse the data. Results The analysis revealed: (1) existence of homogeneity in drug response and resistance development within a patient; (2) tumour shrinkage rate does not relate to rate of resistance development; (3) vemurafenib and dabrafenib, two BRAF inhibitors, have different variability in tumour shrinkage rates. Conclusions Overall these results show how analysis of the dynamics of individual lesions can shed light on the within and between patient differences in tumour shrinkage and resistance rates, which could be used to gain a macroscopic understanding of tumour heterogeneity. Electronic supplementary material The online version of this article (10.1007/s00280-017-3486-3) contains supplementary material, which is available to authorized users.


Introduction
Tumour heterogeneity at the molecular level is known to exist not only between patients but also between lesions within a patient and within an individual lesion [1][2][3]. At the individual lesion level, we can envisage that the molecular heterogeneity is likely to lead to differential cell killing, under a given treatment, within the lesion [4]. The differential killing is likely to vary across lesions within a patient and also across patients. This variability in cell killing could well be visible at the whole tumour level via measurements obtained through routine clinical imaging. Data from clinical trials are likely to be the best source for exploring the variability described as the imaging data collection process is standardised for a large number of patients. This is due to most clinical trials employing the Response Evaluation Criteria In Solid Tumours (RECIST) [5]. This criterion, however, limits the analysis of variability to the between patient and within patient level, as we now explain.
In most cancers, the response of a tumour to treatment is predominantly measured through quantifying images taken of it over time. A standardised methodology to quantify patient response to a treatment that is routinely applied in clinical trials is RECIST. The criterion involves taking information at the individual lesion level and combining it to produce a single value, response category, at each imaging visit for a patient in the following way. The tumours within a patient are first classified as either target or non-target lesions based on whether a lesion is repeatedly measurable or not. Target lesions are then recorded quantitatively by taking the longest diameter of each of them and summing them together to produce the sum of longest diameters (SLD) value at each imaging visit. Non-target lesions are not recorded quantitatively but are recorded qualitatively by assessing whether they have disappeared, still visible/partially shrunk but not grown or experienced unequivocal growth. The information on target lesions, non-target lesions and whether a new lesion has occurred is combined, such that at each visit the patient is placed into one of the following four categories: complete response (CR), partial response (PR), stable disease (SD) or progressive disease (PD). This categorisation scheme is further simplified by combining CR and PR to create the objective response rate (ORR). Once a patient enters the PD category or the patient dies, all imaging stops. This brief introduction to the RECIST criteria clearly highlights that quantitative information at the individual lesion level, through measurement of the target lesions, is available. It is this information that can be leveraged to explore quantitatively the dynamics of response and resistance of tumours at both the between patient and within patient level.
The goal of this study is to explore the variability in the dynamics of the time-series of these target lesions under treatment. This will be done by placing a mathematical model of tumour growth within a statistical framework used routinely for population analysis. The model and statistical framework will allow us to explore the following three biological questions: (1) is there a degree of correlation in the dynamics of tumour size within a patient; (2) if so what is the difference in-between and within patient variability in tumour shrinkage and resistance rates; (3) is there a correlation between tumour shrinkage and resistance rates at either the individual lesion level or patient level. The framework described in this study is applied to three treatments currently used within the metastatic melanoma setting, vemurafenib (BRAF/CRAF inhibitor) [6,7], dabrafenib (BRAF inhibitor) [8] and trametinib (MEK inhibitor) [9]. These compounds were chosen as the pathway under target is the same, and all three are used within the same patient population [10].

Patients
Data from the vemurafenib [11], dabrafenib [12] and trametinib [13] arms of their corresponding phase III studies was collected through clinicalstudydatarequest.com. For full details of the studies and patient demographics, we refer the reader to the previous three articles, which published the results of the phase III studies. Only patients who had a SD, PR or CR response at the first visit were taken forward for the model-based analysis, as our interest is in the response to the treatment followed by tumour resistance.

Derivation of mathematical model of tumour growth
The choice of growth law to be used to analyse the data was based on prior knowledge of our understanding of tumour growth based on empirical observations and biological understanding.
If cells have a cell cycle length t d , then the total number of growing cells will double every t d hours, so their volume will be given by where While this suggests that tumour volume will grow in an exponential or modified-exponential fashion [14], it has often been observed empirically that tumour diameters, as opposed to volumes, appear to grow in a roughly linear fashion. Indeed, this has been known since at least the 1930s. As Mayneord [15] proposed, it was because growth was concentrated in an outer layer of proliferating cells, with cells inside that layer necrotic or quiescent.
Following Mayneord, if we assume that the proliferating layer has thickness d, which is assumed to be small relative to the radius r, and is growing at a rate a, then the volume of the layer is approximately (see Fig. 1). and it is growing at a rate The total volume of the tumour is To translate from cell population growth (with growth rate a) to tumour growth, we, therefore, need just two additional parameters, which are the thickness of the growing layer d, and the initial radius R 0 . The linear growth in diameter translates to cubic, rather than exponential, growth of the tumour volume.
This idea that tumour growth is driven by an outer layer of proliferating cells, surrounding a quiescent or necrotic core, has been featured in recent mathematical models that simulate treatment of tumours by anti-cancer drugs [16]. From a data analysis perspective, it also offers a number of clear advantages, since it allows the use of easily understood statistical techniques (this is not a justification, but it is certainly a convenience). It also requires a minimal number of parameters, which is appropriate for the analysis of clinical studies that are subject to a high degree of noise and are susceptible to over-fitting.
The linear growth equation will of course not be a perfect fit for the growth of all tumours. It assumes that the thickness of the growing layer is small relative to the overall tumour radius (small tumours will see volume grow in a more exponential fashion). Also, it does not account for the saturation observed at larger volumes, therefore, only applies for tumours of intermediate size.
We assume that the linear equation holds for tumour shrinkage as well as tumour growth; this seems justifiable given that, for the tumour to shrink, dead cells need to be removed, which to a first approximation can be modelled as a linear process. Most relevantly for this study, as with any other simple growth law, it does not account for the effects of resistance. As seen next, we have, therefore, modified it in the simplest way possible, by introducing a separate linear growth rate for the resistant phase.

Individual lesion time-series
Individual lesion time-series were modelled using a piecewise linear model, and this was done in two parts. In the first part, we treated each lesion as being independent from each other, i.e. we did not account for which patient the lesions belonged to. This model was represented by the following pair of equations, where subscript j represents each lesion (j = 1,…, m), subscript k represents each time-point (k = 1,…, n), L jk and e jk are the longest diameter and residual error, respectively, for lesion j at time k, BSL j , d j and g j are the initial longest diameter value, decay rate and re-growth rate, respectively, for lesion j. The switching time-point, sp, represents the switch from decay to re-growth. The value of sp was determined in the following way due to identifiability issues with allowing it to vary across lesions. We created a small set of possible sp values, day 63, 116 and 179, by taking the mid-point between on-treatment imaging visit timepoints, on days 42, 84, 147 and 210. For each sp value, we fitted the pair of linear equations and chose the sp value which gave the best fit according to the log-likelihood; higher log-likelihood implies better fit. This approach gives an approximate switching time-point rather than an exact value.
The second part involved accounting for which patient the lesions belonged to. This was done by modifying the above pair of equations by introducing a new level of hierarchy i, which represents each patient (i = 1,…, p), Both models were placed within a mixed effects statistical framework. This approach assumes that the individual effects (random effects component) do not correlate with the independent variables. Within this analysis framework, the parameter BSL was assumed to follow a log-normal distribution, chosen to ensure positivity of lesion size values, and all other parameters assumed to be normally distributed. Furthermore, all the random variables were assumed to be independent. The within and between patient variability in decay and re-growth rates was explored through the distribution of the model parameters, using the coefficient of variation (standard deviation divided by the mean).
All analyses were done in R v 3.0.2 with the nlme package used for the mixed effects analysis.

Patients and data
The imaging characteristics of all the patients used in the analyses here can be seen in Table 1. The table highlights that in terms of treatment response, either via objective response rate (ORR) or % change in the sum of longest diameters (SLD) at week 6 (when the first on-treatment imaging visit occurred), dabrafenib and vemurafenib showed very similar outcomes compared to trametinib. These findings mirror the full original study results. It is also noticeable that the number of patients is larger in the vemurafenib study than the dabrafenib and trametinib studies; again this mirrors the original studies.
The time-series of the individual longest diameters for all lesions across the three studies can be seen in Fig. 2. It shows that the frequency of data collection is consistent over time and that the distribution of initial values is similar across all studies. Figure 3 shows the number of lesions per patient across the studies; which highlights that 80 percent of patients across the studies have more than one target lesion. Overall, the visual analysis of the imaging data suggest that the patients selected for the time-series analysis were well matched across all three studies with respect to imaging data collection.

Individual lesion time-series analysis
The piecewise linear models for the individual lesion timeseries described in the Methods section were fitted to tumour data, and the final models (used throughout the rest of the study) were chosen based on the higher log-likelihood (see also the Supplementary Tables S1, S2 and S3). The fits to the final piecewise linear model for each study, can be seen in Fig. 4. Each point in the plots represents a pair of values, observed and fitted. All the points in each plot lie close to the line of unity which implies that the final model describes the data well. Notably, the final model for each study included information on which patient the lesions belonged to, suggesting there is a degree of correlation in tumour size dynamics under treatment within a patient.
Having established that the extra information on which lesion belongs to which patient is important, we next explore the between and within patient variability of   a full table of model parameter values, see Supplementary information Table S4.) In regard to the rate at which the tumour shrinks, we find that both within and between patient variability (coefficient of variation) are considerably different for each drug. The variability is highest for vemurafenib, followed by trametinib and finally by dabrafenib (for which the variability can be considered quite low). However, for a given drug, no difference in the between and within patient variability was found. Similarly, for the tumour re-growth rate, we find that different inferences can be made for the different drugs. Notably, no variability in the tumour regrowth rate (between and within patient) was observed for vemurafenib (see Supplementary Table S1 for more details). Moreover, no difference similar to the extent seen within the decay rate between dabrafenib and trametinib was found.
Fitting the best model (as determined by the log-likelihood) to clinical data, we obtained that the switching point between tumour decay and re-growth, at 63 days, is the same for all drugs. It must be noted that this does not imply that the switch from the decaying and re-growth phase is not the same across lesions/patients, just that from a population perspective this was an optimal value. The value of 63 days is in-between the first, week 6, and second, week 12, ontreatment imaging visits; see Supplementary Tables S1, S2 and S3 for log-likelihood values at other time-points.
Overall, these results show that the piecewise linear model highlights both qualitative and quantitative differences between these drugs, when comparing both between and within patient variability of tumour size dynamics.
The final question to address in this study is whether there is any correlation between the decay rate and regrowth rate of tumour lesions. For vemurafenib, no distribution was required for the re-growth rate in the final   Fig. 6 we focus on dabrafenib and trametinib (where a distribution on the re-growth rate was required in the final model; see Fig. 5). The first observation to note is that we have a small number of lesions (~ 4%) under treatment that have either a positive decay rate (i.e. they are growing) or a negative re-growth rate (i.e. they are shrinking). These results are possible because we placed no restriction on the sign of the decay or re-growth rate during model fitting. The second observation is that there is no correlation between tumour shrinkage and re-growth rate (coefficient of determination: r 2 < 0.1). These findings suggest that how quickly a tumour shrinks under treatment has no relationship to how quickly resistance within that lesion occurs for any of the drugs.

Discussion
The analysis of tumour heterogeneity within a patient has been predominantly explored at the molecular level using various genetic techniques [17,18]. Measuring the complete heterogeneity of all the lesions within a patient at the genomic level is clearly a difficult task for which no known accurate well-validated method currently exists. However, exploring heterogeneity in tumour response to treatment via measuring the size of individual lesions over time is achievable using routinely collected clinical trial imaging data [19]. Clearly, this does not provide details on the mechanisms of resistance. However, it may provide details on the behaviour of the resistant phenotypes through analysis of re-growth rates. It may also allow us to look at how heterogeneous drug response within a patient relates to drug response across  patients. The purpose of the analysis conducted here was to explore the dynamics of tumour size using a mathematical model derived based on empirical observations and biological understanding. This model was subsequently placed within a mixed effects statistical framework to analyse the variability in dynamics at the between patient and within patient level.
As an exemplar of the approach we applied it to three phase III studies, vemurafenib (BRAF/CRAF inhibitor) [11], dabrafenib (BRAF inhibitor) [12] and trametinib (MEK inhibitor) [13]. All three drugs target the same pathway, mitogen activating protein kinase (MAPK) but in slightly different ways. The three studies were conducted in a similar patient population, BRAF mutant positive metastatic melanoma patients, with RECIST criteria used for imaging data collection. As our goal was to analyse the dynamics of lesions that have quantitative data over time, we restricted our analysis to patients who had more than a minimum of one on-treatment imaging visit, and to their RECIST defined target lesions. As well as the patient population being similar, data collection and initial distribution of target lesion size were also similar across all three studies.
To analyse the target lesion size time-series data, we derived an empirical model based on biological observations and mechanistic insights on tumour growth. The empirical model, which was piecewise linear, was then placed within a mixed effects framework to analyse the data. The main results of the analysis were: 1. Knowing which lesion belongs to which patient improved our ability to describe the data over assuming all lesions are independent of each other.

The within and between patient variability in decay rates
is different for all three drugs. 3. For the re-growth rate, no variability was required for vemurafenib, and the within and between patient variability was similar for trametinib and dabrafenib. 4. Finally, no correlation was found between target lesion shrinkage and resistance growth rate.
The first result suggests that although there is heterogeneity in the dynamics of tumour size over time, there is also some degree of correlation in these dynamics across lesions within a patient. This suggests that there may be a degree of homogeneity across lesions, which is likely the result of the system-level effect that the drug has on the patient (i.e. inactivates the same components of the MAPK pathway across all lesions). This system-level effect seems to be more pronounced in dabrafenib and trametinib (low tumour shrinkage variability) and less pronounced in vemurafenib (high tumour shrinkage variability). These conclusions are relevant in the context of recent studies [20] which showed an improved overall survival in metastatic melanoma with dabrafenib (18.2 months) and trametinib (15.6 months) versus the survival with vemurafenib (13.6 months).
The differences seen across treatments, in terms of shrinkage and resistance dynamics (points 2 and 3), could be attributed to the fact that these drugs have different pharmacological profiles. Indeed, there have been reports highlighting the subtle difference between the two BRAF inhibitors, vemurafenib and dabrafenib, both preclinically and clinically [21][22][23]. Overall, the differences highlight that drugs that appear to give the same study level results, vemurafenib and dabrafenib, can be differentiated at the individual lesion level.
The final result was that there was no correlation between the rate of tumour shrinkage and the growth rate of the resistant clone. This result is quite important as it highlights that the rate at which drug-sensitive cells are killed has no bearing on how quickly a resistant clone will grow. This result could have implications for how these drugs and maybe new treatments are dosed. That is there may be no need to use doses and schedules that aim to eradicate tumour cells quickly. This could lead to treatments being less toxic to a patient and increase the options for combination therapies.
This study is not without its caveats. The model used to analyse the dynamics was unable to accurately estimate the switching point from decay to re-growth due to identifiability issues related to the sparseness of the data. Thus, exploring the variability in this switching point across lesions and patients was not possible and is likely to be an issue with other clinical imaging datasets too given the study design here is typical of the field.
In summary, the mathematical modelling and analysis approach undertaken here highlights how more information can be gained from routinely collected clinical data. We hope this encourages the community to consider analysis at the individual lesion level in addition to the patient level results that are routinely reported.