A Biophysical Systems Approach to Identifying the Pathways of Acute and Chronic Doxorubicin Mitochondrial Cardiotoxicity

The clinical use of the anthracycline doxorubicin is limited by its cardiotoxicity which is associated with mitochondrial dysfunction. Redox cycling, mitochondrial DNA damage and electron transport chain inhibition have been identified as potential mechanisms of toxicity. However, the relative roles of each of these proposed mechanisms are still not fully understood. The purpose of this study is to identify which of these pathways independently or in combination are responsible for doxorubicin toxicity. A state of the art mathematical model of the mitochondria including the citric acid cycle, electron transport chain and ROS production and scavenging systems was extended by incorporating a novel representation for mitochondrial DNA damage and repair. In silico experiments were performed to quantify the contributions of each of the toxicity mechanisms to mitochondrial dysfunction during the acute and chronic stages of toxicity. Simulations predict that redox cycling has a minor role in doxorubicin cardiotoxicity. Electron transport chain inhibition is the main pathway for acute toxicity for supratherapeutic doses, being lethal at mitochondrial concentrations higher than 200μM. Direct mitochondrial DNA damage is the principal pathway of chronic cardiotoxicity for therapeutic doses, leading to a progressive and irreversible long term mitochondrial dysfunction.


S1 ETC Inhibition
To model ETC inhibition by DOX, dose dependent Hill functions were constructed for each of the ETC complexes using IC 50 values reported in the literature [3]. The Hill coeficient was determined by digitizing the figures in [3], fitting a Hill function to the data obtained and rounding to the nearest integer. The obtained value was nH = 3. Section S1.1 presents tests varying this parameter. These functions, depicted in Figure S1, were used to scale the activity of the ETC: 3 (1) Where [DOX] is the DOX concentration in µM . Following the ETC inhibition strategy dicussed in section S1.2, these inhibition functions were introduced in the ETC Equations 31, 32, 53, 72, 73, 78, 79 and 94 of the mitochondrial model listed in section S4.2.

S1.1 Variations in nH
The Hill coeficient, nH, presented in Section S1, was adjusted based on digitized data from the figures of [3] and rounded to the nearest integer. Here, tests were performed to see the effects of variations in this parameter on the results. Figure  S2 presents the acute effects of ETC inhibition assuming different values for nH. In this Figure we can observe that the results present the same qualitative behaviour with a small variation in the concentration in which the mitochondrial function collapses. For the simulations of the chronic effects of DOX, as the concentrations were up to 30µM , which is much lower than the reported IC 50 values for the ETC, the results were unaffected by variations in nH. Therefore, different choices for the value of this parameter do not alter the conclusions of this study.

S1.2 Strategies for ETC inhibition
Two different strategies were tested to model ETC inhibition. Strategy 1 was inherited from the original mitochondrial model adopted [2]. This strategy involves ETC inhibition by using the functions described by equations 1 to 4, to scale specific rate constants in the model. Strategy 2, was to inhibit the ETC by using the functions described by equations 1 to 4, to scale the total activity of each complex of the ETC. Both strategies generated qualitively similar results, as depicted in Figure S3. For the simulations of the chronic effects of DOX, as the concentrations were up to 30µM , which is much lower than the reported IC 50 values for the ETC, the results were unaffected by the choice of ETC inhibition strategy. Therefore, the choice of ETC inhibition strategies do not alter the conclusions of this study. So, all other results in this study were generated with strategy 1, to remain consistent with the adopted mitochondrial model.  Figure S3: Simulated effects of ETC inhibition taking into account different strategies.

S2 Redox Cycling
In our model redox cycling was considered to be proportional to the concentraion of DOX in the mitochondria. This effect was incorporated in the model by increasing the O − 2 production in Complex I, which has been identified as the redox cycling site for DOX [4]: Where a 42 is the rate of O ,− 2 production by Complex I, [DOX] is the DOX concentration in µM and k RC is a redox cycling constant. A complete list with all equations of the ETC model adopted in this work can be found in the end of this supplemetnal material.
The redox cycling constant, k RC = 1/15µM −1 , was adjusted using experimental data from Group E of [5] which is the only group in which measurements were made while the drug was still present in the system. This experiment showed a 7% increase in the O .− 2 concentration in rats 2 hours after treatment with 1mg/kg of DOX , which is the equivalent of a 30µM dose in our model.   Figure S4: Superoxide concentration variation for a single therapeutic DOX dose. The red errorbars are experimental measurements of superoxide concentration for untretead rats and 2h after administration of 1mg/kg of DOX.

S3.1 Protein expression
In our model, the expression of the mtDNA encoded proteins, mt prot , was considered to be proportional to the mtDNA content. The protein densities from the original model, ρ C1 , C3 tot , C4 tot and ρ F 1 , corresponding to Complex I, III, IV and ATP synthase respectively, were scaled using a function to relate mt prot and mtDNA content: The constants, k b = 1.42 and k prot = 0.32, were fitted to experimental data [5] where paired measurements following exposure to DOX, showed a correlation between the mtDNA content and the activity of mtDNA encoded proteins, as shown in Figure S5.  Figure S5: Relationhip between mtDNA content and mtDNA encoded proteins expression in our model [5].

S3.2 Fitting mtDNA damage and repair model
The four parameters in the equation for the mtDNA content time derivative, repeated here for convenience, were adjusted using experimental data and constraint to the baseline conditions of the model.
The function that defines the phase plot of mtDNA in the absence of DOX, depicted in the top panel of Figure 4 is: The first assumption is that the model is in steady state at baseline conditions. To enforce this, we have f (mtDN A b ) = 0, where mtDN A b is the mtDNA content at baseline conditions. As [ • OH] n = 1 at baseline conditions, from Equation 7 we get a hard constrain for one of the parameters: Differentiating equation 7 we obtain: The second assumption is that the mtDNA content is in a stable steady state at baseline. To enforce a this, we set: Combining Equations 9 and 10 we get: Combining Equations 8 and 11 and simplifying we get: By assuming mtDN A b = 0.75 and numerically calculating the hydroxil radical concentration derivative with respect to mtDNA at that point, we get the following soft constraint: After these constraints were applied, a parameter sweep was performed in order to find the parameter set that minimizes the the error between the model predictions and experimental data. The error associated to a parameter set was calculated using the L 2 -norm: Where x 1,exp and x 2,exp are the data points corresponding to Group B and C from [5] respectively and x 3,exp corresponds to Group B from [6]. To calculate the corresponding simulated values, x 1,sim , x 2,sim and x 3,sim , simulations were performed replicating the exprimental protocols of each of the experimental data points. The lower bound for κ is defined in Equation 13 and for α and γ, it was zero by definition. The upper bound for κ was 0.16 as the results presented low sensitivity for values greater than that, with the error monotonically increasing. The upper bounds for α and γ were 1e −10 and 3e −5 respectively, as simulations with values greater than those, generated no results within the reported experimental errorbars. Figure S6 shows the results for all simulated parameter sets, where the colorbar represent the error and the white areas represent parameter sets that generated results outside of the reported experimental errorbars. Table S1 presents the optimal parameter set found which generated an error of e=0.1595. It is possible to observe in Figure S6 that all parameter sets that matched experimental data present a value for γ equal or greater than the value of the parameter set used in the study. Therefore, for all of these parameter sets, direct damage to mtDNA by DOX is the primary toxicity pathway responsible for triggering the vicious cycle that leads to mitochondrial dysfucntion as concluded in our study. The sensitivity of the error with respect to the parameters was calculated using Equation 15.
where S i is the sensitivity of the error with respect to the parameter p i .  Figure S6: L 2 norm errors of the parameter sets tested. White areas correspond to parameter sets which generated results outside of the errorbars in at least one point.

S3.3 Hydroxil radical production
This section describes how the normalized concentration of hydroxil radical, [ • OH] n , in Equation 6 is calculated in the model. The model assumes that hydroxyl radical is produced in a two step iron catalyzed Haber-Weiss reaction. In the first step of the Haber-Weiss reaction, Equation 16, superoxide oxidizes proteins containing F e 3+ in iron-sulfur clusters, inactivating these proteins and releasing F e 2+ [7]. In the second step of the Haber-Weiss reaction, Equation 17, also known as Fenton's reaction, the free F e 2+ then reacts with peroxide, producing hydroxil radical: The model assumes that the concentrations of O 2 , and of the catalyzer F e 3+ are constant, as these concentrations are much greater than the concentration of superoxide. Therefore, the hydroxil radical production can be calculated using the net reaction, without explicitly calculating the concentration of the intermediate reactant F e 2+ : From Eq 18, the rate of hydroxyl radical production, j p • OH is: Hydroxil radical is very reactive and upon formation quickly reacts indiscriminately with substances within a 93Å radius with a half-life of 10 −9 s [8]. Therefore, the model assumes that the concentration of hydroxil radical concentration is in steady state, i.e. the rate of production, j p • OH , is equal to the rate of consumption of hydroxil radical j c • OH : Where [T] is the combined concentraion of all molecules that are targets of hydroxil radical oxidation, which include multiple proteins and lipids. As the concentration of hydroxil radical is much less than the combined concentration of oxidation targets, [ • OH] [T ], we assume that [T] is constant and that the consumption of hydroxil radical is directly proportional to its concentration: Where k 2 = k c · [T ]. Substituting Equation 21 in Equation 19 and normalizing with respect to baseline conditions we get: This concentration, [ • OH], is then used in Equation 6 to calculate the term related to oxidative damage to mtDNA.

S3.4 Iron Chelating therapy
The co-administration of dexrazoxane along with DOX has been shown to prevent the rise in free iron levels observed when administring DOX in isolation [9]. Therefore, in this study, a simplified model for the effect of dexrazoxane co-treatment is adopted by assuming that the concentration of F e 2+ is kept at baseline during the iron chelating treatment. So, during dexrazoxane treatment, from Equation 17 we get: Following the same steps from section S3.3, we get that during dexrazoxane treatment the hydroxil radical concentratil is given by: S4 Model Description

S4.1 Parameters
All simulations in this work consider that the mitochondria are in state 3 respiration with the following input parameters:  Figure S7: Reaction schemes for the ETC model adopted.
a 57 = a * 57 · (QH 2 ) 1/2 · C1 inhib (32) Where a xy is the transition rate from state x to y. Superoxide production by Complex I involves a single electron reaction, in the original mitochondria model this was represented as a two electron reaction. We have updated this reaction, represented in Equation 34, to reflect this.
Flux equations: Where F x is the occupancy probability of state x. Complex II/III: System of differential equations describing complex III: Rate equations for complexes II and III: Reaction 1 -reduction of Q to QH 2 : Reaction 2 -diffusion of QH 2 across the membrane: Reaction 9 -FeS to cyt c 1 e − transfer Reaction 10 -superoxide production from Q .− p k 10 = k010 · Keq10 (83) Reaction 33 -cyt c 1 to cyt c e − transfer v 33 = k 33 · cytc1 red · cytc ox − k −33 · cytc1 ox · cytc red (87)