Systems modeling accurately predicts responses to genotoxic agents and their synergism with BCL-2 inhibitors in triple negative breast cancer cells

Triple negative breast cancer (TNBC) is an aggressive form of breast cancer which accounts for 15–20% of this disease and is currently treated with genotoxic chemotherapy. The BCL2 (B-cell lymphoma 2) family of proteins controls the process of mitochondrial outer membrane permeabilization (MOMP), which is required for the activation of the mitochondrial apoptosis pathway in response to genotoxic agents. We previously developed a deterministic systems model of BCL2 protein interactions, DR_MOMP that calculates the sensitivity of cells to undergo mitochondrial apoptosis. Here we determined whether DR_MOMP predicts responses of TNBC cells to genotoxic agents and the re-sensitization of resistant cells by BCL2 inhibitors. Using absolute protein levels of BAX, BAK, BCL2, BCL(X)L and MCL1 as input for DR_MOMP, we found a strong correlation between model predictions and responses of a panel of TNBC cells to 24 and 48 h cisplatin (R2 = 0.96 and 0.95, respectively) and paclitaxel treatments (R2 = 0.94 and 0.95, respectively). This outperformed single protein correlations (best performer BCL(X)L with R2 of 0.69 and 0.50 for cisplatin and paclitaxel treatments, respectively) and BCL2 proteins ratio (R2 of 0.50 for cisplatin and 0.49 for paclitaxel). Next we performed synergy studies using the BCL2 selective antagonist Venetoclax /ABT199, the BCL(X)L selective antagonist WEHI-539, or the MCL1 selective antagonist A-1210477 in combination with cisplatin. In silico predictions by DR_MOMP revealed substantial differences in treatment responses of BCL(X)L, BCL2 or MCL1 inhibitors combinations with cisplatin that were successfully validated in cell lines. Our findings provide evidence that DR_MOMP predicts responses of TNBC cells to genotoxic therapy, and can aid in the choice of the optimal BCL2 protein antagonist for combination treatments of resistant cells.


Mathematical modeling overview
This supplementary modelling overview was used from Lindner et al (1) and adapted to include the new BCL2 inhibitors section. The signaling pathway of the BCL2 family proteins was modeled by a pseudo-reaction network. Using mass-action kinetics, this network was subsequently translated into a set of ordinary differential equation (ODEs) which describe the changes in concentration of a single protein or a single protein complex at a given time point t (denoted as 'reaction rate', or 'rate', hence further). ODEs were parameterized using cell specific concentrations and kinetic constants as input. Subsequently, ODEs were solved using MATLAB 7.3 (MathWorks, USA, R2007b, 7.5.0.342) and its function ode15s.
Lucantoni et al 2 We translated the modeled protein-protein interaction from a general reaction notation to ODEs by a detailed guideline listed in Supplementary Table 1. The guideline is grouped into different reaction classes. In a nutshell, the degradation rate of a protein was calculated by multiplying the actual protein concentration by the degradation constant k deg and this rate was subtracted from the actual protein concentration (Supplementary Table 1a). Protein turnover was realized by adding a production k prod rate to the degradation rate (Supplementary Table 1b). Proteins that sustain turnover were modeled to reach equilibrium concentrations that were given by the interplay of protein production and protein degradation. As a consequence the equilibrium that was attained was given by the ratio of production rate and degradation constant (k prod / k deg ).  Table 1c). The change in concentration due to a protein-protein interaction, such as binding, was modeled by multiplying the product of the concentration of all interacting proteins by the reaction specific constant k (Supplementary Table 1d). This was irreversible unless a backward equation was modeled with a separate reaction specific constant (Supplementary Table 1e). The ratio of the backward constant k backward and the first / forward constant k forward is defined as dissociation constant K D . A Low value of K D indicates a high affinity of the interacting proteins. Vice versa a high value indicates a low affinity.

Modeling of BCL2 protein reactions
Proteins of the BCL2 protein family were modeled to be located and to interact at the mitochondrial outer membrane (2,3). As only exception cytosolic BAX was modeled to have a cytosolic fraction BAX c (4). 3   As the first type of protein interaction, protein degradation was modeled for all proteins and   protein complexes except for BAK, BAX and VDAC2 and their complexes (Supplementary   Table 2). Degradation constants k deg were calculated with the equation k deg = ln(2)/60*t 1/2 from half-life times t 1/2 . We obtained the protein half-life time of the proteins BCL2, BIM, MCL1, NOXA, PUMA and tBID from literature (5)(6)(7)(8)(9)(10)(11). We considered for BCL(X)L the same half-life time as for BCL2 (300 minutes) due to the similar amino acid sequence of BCL2 and BCL(X)L (EMBOSS Needle alignment, similarity 53%). Because MCL1 is subject to rapid degradation (5,7), the degradation constant for MCL1 was modeled to drop from 45 minutes to 17 minutes upon induction of BIM, PUMA and NOXA stress at time point t = 0h after stress induction. We  Table 3). For the binding of BIM to BCL2, BCL(X)L, and MCL1, the dissociation constants K D (K D is defined as k backward / k forward ) as well as the constants k forward and k backward were taken from (14). For the binding of PUMA, NOXA and tBID to BCL2, BCL(X)L, and MCL1, K D values were taken from literature (4,(15)(16)(17)(18)(19)(20)(21). The backward kinetic constants k backward for the reversible binding of BCL2, BCL(X)L, and MCL1 with each BAK, BAX, Lucantoni et al 4 NOXA, PUMA and tBID were modeled to be the same as those for BIM. The forward dissociation constant was calculated by k forward = k backward / K D .

Lucantoni et al
Inactive BAX c was modeled to be only present in the cytosol. Hence, no interaction of inactive BAX c with anti-apoptotic BCL2 family proteins (22) was modeled. BAX c activation by BH3-only proteins BIM, PUMA (and tBID) was modeled in a two-step procedure. Therefore BIM, PUMA (and tBID) were first modeled to reversibly bind with an fitted dissociation constant K D of 100 nM that is consistence with Dai et al. (23) and forward kinetic constant k forward of 2.57x10 -6 nM -1 s -1 . The backward constant was calculated by k backward = K D * k forward . Formed activator~BAX c hetero-dimers were modeled to irreversibly activate BAX c * followed by an instant release of BAX c * with an estimated half-life time of 1/10 minute (Supplementary Table 4a and b). Upon its activation, BAX c * was modeled to translocate into the mitochondrial outer membrane (MOM) with a fast single cell kinetic that was determined by our lab previously (24) (Supplementary Table 4c).
In the model, BAK activation was modeled with the same reactions and kinetics (Supplementary Table 4a and b) as BAX. However, these interactions were considered to take place at the mitochondrial membrane and therefore no translocation was modeled (25). Unlike BAX, BAK was modeled to interact with the Voltage-dependent anion-selective channel protein 2 (VDAC2) (26,27). We estimated that VDAC2 was expressed with the same concentration as BAK in the respective cell. Only inactive BAK was modeled to bind to VDAC2. The VDAC2~BAK dimer was modeled to dissociate with an estimated backward constant k backward = 2.31x10 -3 s -1 and a dissociation constant of 1,000 nM. The forward constant k forward was calculated to be 2.31x10 -6 nM -1 s -1 (Supplementary Table 4d).
Lucantoni et al 5 Active BAK * and BAX * was modeled to homo-oligomerize up to dodecamers (Supplementary Table 5). For homo-oligomerization, a K D of 15 nM and a decay time of the pores of 1 hour were estimated. Hexamers or larger homo-oligomers were modeled as mitochondrial pores. Once 10% of total effectors formed pores, mitochondrial outer membrane permeabilization (MOMP) was considered to be induced, in analogy with previous studies from our group (24).
To study the effect of apoptosis sensitizers, further ODEs that characterized the interaction of the BH3-mimetics ABT199, WEHI-539 and A-1210477 with the BCL2 family proteins (Supplementary Table 6) were included. Each of these drugs was modeled to bind to the anti-apoptotic proteins BCL2, BCL(X)L and MCL1. Specific dissociation constants K D were taken from the literature (28)(29)(30)(31), the backward constants k backward were fixed to reasonable values and the forward constant k forward resulted from k backward and K D . The degradation constants k deg of the antagonists and their hetero-dimers with BCL2, BCL(X)L and MCL1 were set equivalent to the degradation constants of the hetero-dimers of anti-apoptotic BCL2 proteins and the BH3-only proteins, respectively. These constants remained unchanged, as described in Lindner el al. (1) .

Model Input
We used two major inputs for the model. First we used absolute protein concentrations of BAK, BAX, BCL2, BCL(X)L and MCL1. Secondly, we modeled cell stress by protein production of BH3-only proteins.
We determined the concentration of BAK, BAX, BCL2, BCL(X)L and MCL1 by quantitative Western Blotting (see Figure 2 and Table 1). Since BCL2, BCL(X)L and MCL1 sustain protein degradation with half-life times of less than or equal to 5 hours (5-11), we considered a turnover for these proteins. The cell and protein specific production rate k prod were calculated by multiplying the quantified (target) protein concentration by the respective degradation constant  Table 2). We modeled stable levels of BAK and BAX (protein half-life time > 22 hours (12,13)) and therefore modeled neither protein production nor protein degradation for those proteins.
Stress induction was modeled by protein production of BIM/PUMA/NOXA (32). Unless noted otherwise, the duration of BH3-only production rate was 12 hours, between time points t = 0h and t = 12h after stress induction, and kept on a constant level. The BH3-only stress dose η is defined as the production rate per hour multiplied by the duration. Where indicated, the BH3-only stress dose η required for MOMP of these proteins was calculated by using the iterative procedure as described below. Like the BH3-only stress dose η, administration of ABT199, WEHI-539 and A-1210477 was modeled with a continuous increase to a specific total concentration of the respective drug over 12 hours (33,34), between time point t = 0h and t = 12h after stress induction.
For each cell, a steady state of the protein concentration was calculated in the absence of stress.
The resulting steady state was used as the initial state for the subsequent calculation in the presence of stress.

Determining of the model predicted BH3-only stress dose η that is required for MOMP
The BH3-only protein stress dose η (BIM/PUMA/NOXA) that is required for MOMP was determined by an iterative approximation algorithm for each cell. MOMP was considered to occur once 10% of total effectors were bound to pores (35). An initial protein stress rate (protein production per time, ν 0 ) of 2 μM/h (correspond to a 24 μM dose over 12 hours) was used. The initial step size for the iteration Δν 0 was set to the half value of the initial protein stress rate ν 0 . It was then determined whether or not the model predicted MOMP in that particular cell, under the Lucantoni et al 7 initially used stress. If MOMP was not predicted, the used protein stress rate was increased by the step size Δν i (ν i+1 = ν i + Δν i ). When the protein stress was sufficient to cause MOMP, the used protein stress was decreased by the step size Δν 0 (ν i+1 = ν i -Δν i ). Subsequently, a new step size was set to the half of the actual value of Δν i+1 = Δν i / 2. Calculations were repeated until the step size was smaller than 10 -6 μM/h (correspond to a minimal dose of 1.2x10 -5 μM over 12 hours). Finally, the BH3-only protein stress dose η required for MOMP was calculated by multiplying the determined protein stress ν by the duration for which the stress was modeled to be present (12 hours).
Lucantoni et al 8 Supplementary Table Legends   Supplementary Table 1 (12,13)) and therefore modeled protein neither production nor degradation for

Supplementary Table 3: Pseudo-reactions and kinetics for inhibition of BH3-only proteins and effectors BAK and BAX by anti-apoptotic proteins
Pseudo reactions and kinetic constants for binding of (a) BCL2, (b) BCL(X)L, (c) MCL1 to BH3-only proteins, and BAK and BAX are given. Dissociation constants K D were taken from literature as indicated (14). Whenever no binding was reported, a K D of 10,000 nM was modeled. Wherever indicated, the backward constant k backward was taken from (14) and otherwise fixed as stated in the text. No K D for the binding of BCL2 to tBID was reported in literature. We modeled a K D of 4 nM with respect to the assumption that the ratio of the K D of the binding of BCL2 to BID (66 nM (18)) and BCL2 to tBID is about the same as the ratio of the K D binding of BCL(X)L to BID (448 nM (34)) and BCL(X)L to tBID (27.2 nM (21)). Forward constants were obtained from the backward constants k backward and from dissociation constants K D , by

Supplementary Table 4: BAK and BAX activation and BAK inhibition by VDAC2
Activation of BAK and BAX by BIM, PUMA, tBID and by the tBID chimeras tBID BAK and tBID BAX is given. was modeled to translocate with a half-life time 1/10 minute (k = 1.16x10 -1 1s -1 ) to the mitochondrial membrane. (d) Inactive BAK was modeled to be able to bind to VDAC2.

Supplementary Table 5: Effector homo-oligomerization
Activated BAK and BAX were modeled to homo-oligomerize in the mitochondrial outer membrane. Effectors were able to build homo-dimers which further were modeled to build larger homo-oligomers. Independent of the complex size, the same dissociation constants and kinetic constants were used. Homo-oligomers larger or equal than hexamers were considered as pores. Only homo-oligomers smaller than or equal to dodecamers were taken into account. Size of effector oligomers are indicated by superscripts.

1210477.
Interaction kinetics of the apoptosis sensitizers ABT199, WEHi-539 and A-1210477 with BCL2 proteins were included into the model.

Supplementary Tables
Supplementary Table 1 a protein and protein complex degradation