Quantifying ERK activity in response to inhibition of the BRAFV600E-MEK-ERK cascade using mathematical modelling

Background Simultaneous inhibition of multiple components of the BRAF-MEK-ERK cascade (vertical inhibition) has become a standard of care for treating BRAF-mutant melanoma. However, the molecular mechanism of how vertical inhibition synergistically suppresses intracellular ERK activity, and consequently cell proliferation, are yet to be fully elucidated. Methods We develop a mechanistic mathematical model that describes how the mutant BRAF inhibitor, dabrafenib, and the MEK inhibitor, trametinib, affect BRAFV600E-MEK-ERK signalling. The model is based on a system of chemical reactions that describes cascade signalling dynamics. Using mass action kinetics, the chemical reactions are re-expressed as ordinary differential equations that are parameterised by in vitro data and solved numerically to obtain the temporal evolution of cascade component concentrations. Results The model provides a quantitative method to compute how dabrafenib and trametinib can be used in combination to synergistically inhibit ERK activity in BRAFV600E-mutant melanoma cells. The model elucidates molecular mechanisms of vertical inhibition of the BRAFV600E-MEK-ERK cascade and delineates how elevated BRAF concentrations generate drug resistance to dabrafenib and trametinib. The computational simulations further suggest that elevated ATP levels could be a factor in drug resistance to dabrafenib. Conclusions The model can be used to systematically motivate which dabrafenib–trametinib dose combinations, for treating BRAFV600E-mutated melanoma, warrant experimental investigation.

SM1 Formulating the system of differential algebraic equations SM1. 1 The system of reactions (R.1R.36) The system of chemical reactions (R.1R.36), as listed below, describes signalling dynamics in the intracellular BRAFV600EMEKERK cascade, when a cell is subjected to the BRAFinhibitor dabrafenib and the MEKinhibitor trametinib. The cascade and drug actions are schematically illustrated in Figures 1 and 2 in the main manuscript, respectively, and (R.1R.36) are listed in Figure 3. The dotnotation (·) is here used to denote complexes of two or more molecles. Reactions involving dabrafenib and trametinib are respectively marked by single (*) and double (**) stars, and can be omitted if one wishes to study the system in the absence of these drugs.

SM1.2 The system of ordinary differential equations (O.1O.36)
Using the law of mass action, the system of reactions (R.1R.36) can be formulated in terms of a system of ordinary differential equations (O.1O.36). We here let y 36x1 denote a column vector with 36 elements, where each element corresponds to a signalling molecule concentration so that,  Using this notation, the system of ODEs (O.1O.36) can be formulated as, + y(8)k 12 + y(10)k 12 + y(8)k 12 + y(10)k 12 + y(1)y(15)(−a 4 ) + d 4 y(16), Note that, we assume a cellular ADPtoATP recycling to continually occur and thus the [ADP ] variable is not included in the [AT P ] conservation law. This ATPtoADP recycling is not explicitly included in the model. In C.1C.7, the total concentrations, denoted by a totsubscript, are equal to the initial concentrations described in Table ST2 (Section SM3.2) so that,

SM2 Code access and instructions
The computational MATLAB [1] code includes functions to convert the system of reactions (R.1 R.36) into a system of differential algebraic equations, i.e., (O.1O.36) with codegenerated substi tutions (C.1C.7), using the law of mass action. The MATLAB function ode15s is then used to solve for y 36x1 (i.e., the signalling molecule concentration vector) as a function of system concentrations and time.

SM2.1 Code access GitHub repository
The MATLAB code, developed for this project, is available on the code hosting platform GitHub, in the project repository (https://github.com/SJHamis/MAPKcascades).

SM2.2 Instructions on how to run the code
The model outputs/result plots presented in the main manuscript can be generated by running the five mainfiles, with file names main_[insert option].m. The code includes three subdirectories: (1) auxiliary_files_plots, which includes files needed to generate data for specific time points, (2) auxiliary_files_model_setup, which includes files needed to set up the DAE MATLAB func tion file (mapk_cascade_DAE.m), (3) model_details, in which the system of reactions (R.1R.36) and the rate constants are manually set.

SM2.3 Instructions on how to modify the code
Instructions on how to edit the cascade structure, rate constants and initial conditions are provided below. For more detailed information, see inline comments in the MATLAB code. Figure SF1: Code directories.

SM2.3.1 Modifying the rate constants
The rate constants can be manually edited in the file model_details/set_reaction_constants.m. By first editing the constants and thereafter running the file auxiliary_files_model_setup/main_setup_convert_SoRs_to_SoODEs.m, the MATLAB DAE function file auxiliary_files_model_setup/mapk_cascade_DAE.m will be updated to include the modified rate constants.

SM2.3.2 Modifying the chemical reactions/cascade structure
The system of reactions can be manually edited in the file model_details/set_system_of_reactions.m. By first editing the system of reactions, by for example by including more terms, and thereafter running the file auxiliary_files_model_setup/main_setup_convert_SoRs_to_SoODEs.m, the function files auxilary_files_model_setup/mapk_cascade_DAE.m and auxilary_files_model_setup/get_conslaw_position.m will be updated to include the new model equations (ODEs and conservation laws).

SM2.3.3 Modifying the initial conditions
The model initial condition (i.e., the initial cascade component concentrations) can be modified in the main files, main_[insert option].m.    [2] listed MAPKKK=3nM, and we use this as our baseline value. However, in order to study amplified BRAF levels as a mode of drug resistance, we explore initial BRAF concentrations within the range 310 nM. Similarly, initial ATP concentrations are explored within the range 15mM, but the baseline value is 1mM.

SM4 Data and derivation of model parameters
All model parameter values for the forward rate constants a 1 , a 2 , ..., a 8 , the reverse rate constants d 1 , d 2 , ..., d 8 and the catalytic rate constant k 1,2 , k 3 , k 5,6 , k 7 are obtained from data that is available in the literature. Note that a i = a 1 ∀i = 2, 3, ...8, as is explained in the main manuscript (in Methods Model parameters).

SM4.1 Rate constants for BRAFATP and BRAFMEK reactions
VanScyoc et al. [3] provide kinetic data for the phosphorylation of the substrate MEK by wild type BRAF. We use their data to approximate the parameter values for a 1 , d 2 and k 1,2 in our model, i.e. rate constants for the BRAFATP and BRAFMEK reactions. We also compute the reverse BRAF ATP rate constant for wild type BRAF, here denoted by d 1,W T . We then use d 1,W T to find a 1 , as the catalytic rate constant k 1,2 and the forward rate constant a 1 , can be assumed to be the same for both V600Emutant and wild type BRAF interacting with ATP. In order to set d 1 , we thereafter use data for BRAF V600E ATP interactions provided by Hatzivassiliou et al. [4]. the data for the Michaelis constant to solve for the parameter value of a 1 ,

BRAFATP reactions
(P.3) If we now plug this value of a 1 into the expression for d 1,W T in (P.1) we obtain, (P. 4) In order to now find d 1 , the reverse rate constant for V600E mutant BRAF, we use the ATP BRAF V600E K m value reported by Hatzivassiliou et al. which is 65 µM. Using this, we can obtain d 1 via Note that the catalytic rate constant for BRAFMEK, and ATP, reactions was obtained in (P.2).

SM4.2 Rate constants for MEKATP and MEKERK reactions
Mansour et al. [9] provide steadystate rate constants for MEK reactions. We use data from their study to set model parameter values describing MEKATP and MEKERK reactions. Mansour et al. further report that K m,AT P = 3.5 ± 0.8µM, which corresponds to K m5 in our model so that,

SM4.4 Rate constants for drug reactions
BRAFDBF interactions (approximating d 4 ): Data pertaining to DBF reactions is gathered from the US Food and Drug Administration [5], according to which the in vitro IC50 value is 0.65nM.
Rheault et al. [6] similarly report this IC50 value to be 0.7nM. We use this IC50 value to estimate the inhibitory constant K i via the ChengPrusoff approach so that, where the experimental ligand (ATP) concentration is 10µM [7], and the ATP K d value has been approximated as K m = 65µM, as previously used. We can now find d 4 via, Koelblinger et al. [8] report the reverse rate constant for dabrafenibBRAF interactions to be 9.6 · 10 −5 s −1 , which is in good agreement with our calculated value.