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

Mitogen-activated protein kinase (MAPK) pathways are crucial for cellular processes like proliferation, differentiation, and apoptosis. The BRAF-MEK-ERK cascade is a key MAPK pathway, and mutations in this pathway, particularly the BRAFV600E mutation, are significant drivers of melanoma development and growth. BRAFV600E leads to constitutive activation of the cascade, resulting in uncontrolled cell proliferation. Small-molecule inhibitors targeting BRAF and MEK have been developed to suppress ERK activity and, consequently, cell proliferation. 'Vertical inhibition,' a treatment approach involving simultaneous targeting of multiple tiers in the cascade (e.g., BRAF and MEK), is a standard of care for BRAF-mutant melanoma. While effective, this approach still faces challenges due to the rapid development of drug resistance. Several mechanisms contribute to this resistance, including BRAFV600E amplification and increased ATP levels. This study aims to develop a mechanistic mathematical model to quantitatively understand the synergistic effects of vertical inhibition and the factors leading to drug resistance, ultimately guiding the choice of effective drug combinations and doses.

Previous research established the importance of the BRAF-MEK-ERK pathway in melanoma and the role of BRAFV600E mutations in driving tumor growth. The development of BRAF and MEK inhibitors, such as dabrafenib and trametinib, has revolutionized melanoma treatment, although monotherapies often lead to drug resistance within 6-8 months. Combination therapies involving both BRAF and MEK inhibitors have shown promise in delaying resistance. Studies have identified several resistance mechanisms, including BRAFV600E amplification, which results in sustained ERK activation even in the presence of inhibitors. Existing models like Huang and Ferrell's general MAPK cascade model provide a foundation, but lack the specificity for the BRAFV600E-MEK-ERK cascade and the detailed mechanism of action of specific inhibitors like dabrafenib and trametinib. The current study seeks to address these gaps.

The authors developed a mechanistic mathematical model of the BRAF-MEK-ERK cascade in BRAFV600E-mutant melanoma cells. This model extends Huang and Ferrell's general MAPK cascade model by incorporating specific details of the BRAF-MEK-ERK pathway and explicitly including the actions of the BRAF inhibitor dabrafenib and the MEK inhibitor trametinib. The model accounts for ATP-dependent substrate phosphorylation. The model is represented by a system of 36 chemical reactions, which are then translated into a system of ordinary differential equations (ODEs) using mass action kinetics. The ODEs are solved numerically using MATLAB, incorporating conservation laws to ensure the total concentrations of the kinases, phosphatases, and drugs remain constant. Model parameters, including kinetic rate constants and initial concentrations, were obtained from the literature. The model's output focuses on the time-dependent concentration of doubly phosphorylated ERK (ppERK), which represents activated ERK. Simulations were performed for various concentrations of dabrafenib, trametinib, BRAF, and ATP to assess the effects of monotherapies and combination therapies on ERK activity.

The model successfully quantifies ERK activity under various treatment conditions. Simulations of dabrafenib monotherapy demonstrated a dynamic response highly dependent on BRAF and ATP concentrations. Higher BRAF and ATP levels resulted in reduced sensitivity to dabrafenib. Trametinib monotherapy simulations showed distinct steady-state ERK levels dependent on the drug dose, with minimal sensitivity to ATP variations within the tested range. Combination therapy simulations using dabrafenib and trametinib revealed synergistic effects in inhibiting ERK activity at specific drug concentrations. However, increased BRAF and ATP concentrations significantly reduced the efficacy of combination therapies, indicating a mechanism of drug resistance. Heatmaps visualized the synergistic and antagonistic effects of different drug combinations based on activated ERK levels. The model predicts that higher BRAF and ATP levels are associated with resistance to the combined treatment. A comparison with clinical data from a study by Falchook et al. showed that the model's predictions of ERK activation closely matched patient observations at a BRAF concentration of 3 nM, but deviated significantly when BRAF concentration increased to 10nM, reflecting a drug-resistant state.

The model successfully predicts the synergistic effect of combined BRAF and MEK inhibition and highlights the crucial role of BRAF and ATP levels in mediating drug resistance. The findings align with experimental observations showing that increased BRAFV600E concentrations lead to drug resistance. The model's prediction of ATP's influence on drug sensitivity underscores the importance of intratumoral heterogeneity in treatment response. The close agreement between the model's predictions and clinical data at lower BRAF concentrations validates the model's accuracy. However, the deviation at higher BRAF concentrations suggests the model might need further refinement to fully capture complex interactions at higher concentrations or the interplay between BRAF and ATP concentration. The model's ability to quantify these relationships provides valuable insights for optimizing combination therapies.

This study presents a novel mathematical model to quantify the effect of combined BRAF and MEK inhibition on ERK activity in BRAFV600E-mutant melanoma. The model successfully captures the synergistic effect of the combination therapy and identifies BRAF and ATP levels as key factors in drug resistance. Future work could extend the model to include an ERK inhibitor to simulate three-drug vertical inhibition or incorporate the model into a multi-cellular agent-based model to simulate tumor growth and evolution under various treatment strategies, incorporating spatiotemporal variations in cellular parameters. Further experimental validation of the model's parameters could improve its predictive power.

The model relies on parameters obtained from the literature and might not perfectly capture the complexity of intracellular interactions. The model's accuracy may be limited by the assumptions made in simplifying the system of reactions and by the availability of detailed kinetic data. Further experimental validation with broader parameter ranges is necessary. The model currently lacks explicit representation of the cell cycle, and the relationship between ERK activity and cell proliferation is implied rather than explicitly modeled. Future development may benefit from incorporating cell cycle dynamics and more detailed modelling of the interactions between ERK and cell cycle regulators.