A Mathematical Model for DNA Damage and Repair

Abstract

In cells, DNA repair has to keep up with DNA damage to maintain the integrity of the genome and prevent mutagenesis and carcinogenesis. While the importance of both DNA damage and repair is clear, the impact of imbalances between both processes has not been studied. In this paper, we created a combined mathematical model for the formation of DNA adducts from oxidative estrogen metabolism followed by base excision repair (BER) of these adducts. The model encompasses a set of differential equations representing the sequence of enzymatic reactions in both damage and repair pathways. By combining both pathways, we can simulate the overall process by starting from a given time-dependent concentration of 17𝛽-estradiol (E₂) and 2^′-deoxyguanosine, determine the extent of adduct formation and the correction by BER required to preserve the integrity of DNA. The model allows us to examine the effect of phenotypic and genotypic factors such as different concentrations of estrogen and variant enzyme haplotypes on the formation and repair of DNA adducts.

1. Introduction

Estrogens are known carcinogens recognized as prime risk factor for the development of breast cancer [1]. The principal enzyme controlling estrogen metabolism in breast tissue is cytochrome P450 1B1 (CYP1B1), which sequentially oxidizes E₂ to 4-OHE₂ and the estrogen quinone (E₂-3,4-Q, Figure 1). The highly reactive quinone nonenzymatically attacks DNA and forms covalent adducts with bases, such as 4-OHE₂-N7-Gua (GUA_adduct), which undergoes depurination by hydrolysis of the N-glycosylic bond between base and sugar, leaving an apurinic (AP) site in the double-stranded DNA. AP sites are repaired by the base excision repair (BER) pathway in a sequence of enzymatic reactions. AP endonuclease APE1 cleaves 5^′ to the AP site to generate a nick with 3^′-OH and 5^′-deoxyribose phosphate termini. DNA polymerase (Pol) catalyzes the release of the deoxyribose phosphate residue and DNA synthesis to fill the gap, which is then sealed by DNA ligase (Figure 1). This reaction sequence is carried out in one of two pathways, “short-patch” and “long-patch”, which differ in the number of inserted nucleotides and the protein components.

Figure 1 

Overview of estrogen metabolism-induced DNA adduct formation and subsequent DNA repair by BER. The top sequence of reactions shows the CYP1B1-mediated oxidative estrogen metabolism resulting in estrogen-deoxyribonucleoside adduct formation. The bottom sequence of reactions shows the BER pathway with its two branches, the predominant short-patch pathway (left) and the alternate long-patch pathway (right). Not shown is the formation of the catechol estrogen 2-OHE2, which yields stable adducts, and the redox cycling of catechol estrogens and estrogen quinones, which generates reactive oxygen species capable of producing oxidative DNA adducts.

In this paper, we present a mathematical model for the formation of DNA adducts from oxidative estrogen metabolism followed by BER of these adducts. The modeling is initially done in two stages: (1) formation of estrogen-deoxyribonucleoside adducts, and (2) repair of adducts in the BER. By combining both stages, we can simulate the overall process by starting from a given time-dependent concentration of E₂ and 2^′-deoxyguanosine, determine the extent of adduct formation and the correction by BER required to preserve the integrity of DNA.

2. Mathematical Model

The mathematical model is constructed in stages. Stage 1 deals with the formation of DNA adducts from estrogen metabolism and Stage 2 encompasses the repair of adducts. In the differential equations that will constitute the model, E₂ denotes 17β-estradiol, 4-OHE₂ the catechol estrogen, E₂-3,4-Q the estrogen quinone, GUA the deoxyribonucleoside 2^′-deoxyguanosine, and GUA_adduct the estrogen adduct 4-OHE₂-N7-Gua over time 𝑡. We assume that the starting points for the dynamic model are the initial concentration of GUA(0)=GUA₀ and a prescribed time series is given for E₂.

The first stage of the model is based on the first three steps of the reaction sequence shown in Figure 1. As chemical reactions, these steps areE₂+E_CYP1B1⇄4-OHE₂+E_CYP1B1,4-OHE₂→E₂-3,4-Q,GUA+E₂-3,4-Q→GUA_adduct. (1) In the first reaction, the production of 4-OHE₂ from E₂ is an enzymatic reaction mediated by CYP1B1 that follows Michaelis-Menten kinetics with parameters 𝐾_𝑚1 and 𝑘_cat1 [3]. In the second reaction, E₂-3,4-Q is produced from 4-OHE₂ according to Hill-type kinetics. It has been shown that E₂-3,4-Q is highly reactive and will spontaneously form adducts with DNA without enzyme involvement [2]. Thus, in the third reaction, the formation of adducts is proportional to the levels of the quinone and the deoxyribonucleoside with a rate constant denoted by 𝑘_Q. Estrogen-DNA adducts have been detected in normal and malignant human breast tissues [4, 5], and we have recently provided direct experimental evidence that oxidative metabolism of the parent hormone E₂ leads to the formation of 4-OHE₂ and deoxyribonucleoside adducts, such as 4-OHE₂-N7-Gua [2]. Using a compartmental analysis we created a mathematical model for reaction sequence (1). Figure 2 shows superimposed the experimental data (dots) and the model simulations (curves) for E₂, 4-OHE₂, and GUA_adduct.

(a)

(b)

(c)

Figure 2 

Comparison of stage 1 mathematical model (DNA damage) with experimental data. The metabolism of E2 (a), 4-OHE2 (b), and GUAadduct (c) is shown as a function of time. In each graph the dots represent the experimental data [2], and the curves are derived from the mathematical model.

The second stage of the model is based on the BER pathway (Figure 1) and uses a compartmental analysis with Michaelis-Menten kinetics to simulate the enzymatic reactions in the pathway [6]. The enzymes include the endonucleases APE1 and FEN1, DNA polymerases 𝛽 and 𝛿 (Pol𝛽, Pol𝛿), and ligases 1 and 3 (Lig1, Lig3). The dynamic process is modeled by a system of nonlinear differential equations. 𝑆(𝑡) denotes the concentration of DNA adducts. We define the source and removal function for the repair cycle functions, 𝑆(𝑡) and 𝑅(𝑡), as 𝑆(𝑡)=𝑘_𝑠GUA_adduct and 𝑅(𝑡)=𝑘_𝑟(E₂-3,4-Q)𝑦₉. Here, 𝑘_𝑠 is a modeling parameter expressed in [1/min] which reflects the rate at which GUA_adduct enters the repair cycle. The parameter 𝑘_𝑟 is the rate at which 𝑦₉(𝑡) interacts with E₂-3, 4-Q to form GUA_adduct. The two parameters, 𝑘_𝑠 and 𝑘_𝑟, can be used to control the rates of entering and exiting of adducts and repaired lesions, respectively. For the simulations, we have chosen, 𝑘_𝑟=𝑘_𝑠=1,𝑑𝑦₄𝑘𝑑𝑡=𝑆(𝑡)−^cat5E_APE1𝑦₄𝐾_𝑚5+𝑦₂,𝑑𝑦₅=𝑘𝑑𝑡^cat5E_APE1𝑦₄𝐾_𝑚5+𝑦₄−𝑘_cat6E_Pol𝛽𝑦₅𝐾_𝑚6+𝑦₅−𝑘_cat9E_Pol𝛿𝑦₅𝐾_𝑚9+𝑦₅,𝑑𝑦₆=𝑘𝑑𝑡^cat7E_Pol𝛽𝑦₇𝐾_𝑚7+𝑦₇−𝑘_cat6E_Pol𝛽𝑦₆𝐾_𝑚6+𝑦₆,𝑑𝑦₇=𝑘𝑑𝑡^cat6E_Pol𝛽𝑦₅𝐾_𝑚6+𝑦₅−𝑘_cat7E_Pol𝛽𝑦₇𝐾_𝑚7+𝑦₇,𝑑𝑦₈=𝑘𝑑𝑡^cat6E_Pol𝛽𝑦₆𝐾_𝑚6+𝑦₆−𝑘_cat11E_Lig1𝑦₈𝐾_𝑚11+𝑦₈−𝑘_cat8E_Lig3𝑦₈𝐾_𝑚8+𝑦₈,𝑑𝑦₉=𝑘𝑑𝑡^cat11E_Lig1𝑦₁₁𝐾_𝑚11+𝑦₁₁+𝑘_cat8E_Lig3𝑦₁₁𝐾_𝑚8+𝑦₁₁+𝑘_cat11E_Lig1𝑦₈𝐾_𝑚11+𝑦₈+𝑘_cat8E_Lig3𝑦₈𝐾_𝑚8+𝑦₈−𝑅(𝑡),𝑑𝑦₁₀=𝑘𝑑𝑡^cat9E_Pol𝛿𝑦₅𝐾_𝑚9+𝑦₅−𝑘_cat10E_FEN1𝑦₁₀𝐾_𝑚10+𝑦₁₀,𝑑𝑦₁₁=𝑘𝑑𝑡^cat10E_FEN1𝑦₁₀𝐾_𝑚10+𝑦₁₀−𝑘_cat11E_Lig1𝑦₁₁𝐾_𝑚11+𝑦₁₁−𝑘_cat8E_Lig3𝑦₁₁𝐾_𝑚8+𝑦₁₁. (2) There are several Michaelis-Menten parameters, 𝑘_cat𝑖 and 𝐾_𝑚𝑖, for which experimental values have been estimated or measured (Table 1) [6]. The enzyme levels are denoted by E_APE1, E_FEN1, E_Pol𝛽, E_Pol𝛿, E_Lig1, and E_Lig3. The variable, 𝑦₉(𝑡), represents the repaired GUA_adduct, and if we set its initial condition in the system of differential equations as 𝑦₉(0)=GUA₀, it has the same meaning as GUA(𝑡). AP endonucleases, such as AP endonuclease 1 (APE1), are quite abundant in most cells, ranging from 200,000 to 7,000,000 APE1 molecules per cell in fibroblasts and HeLa cells, respectively [7, 8]. Thus, the minimal concentration of APE1 in nuclear extract can be assumed to be 2000nM [6]. The equivalent estimated concentrations of Pol𝛽, Pol𝛿, FEN1, Lig1, and Lig3 are 419, 600, 450, 254, and 254 nM, respectively [6, 9, 10].

Table 1 

Parameters for System of Differential Equations [6].

It is estimated that AP sites are also generated spontaneously at an estimated rate of 2,000 to 10,000 per cell per day [11, 12]. We used (2) to determine the time necessary to repair these AP sites. In Figure 3 we demonstrate that the repair of lesions is a nonlinear process. For small values of the ratio (APlesions/10⁵Nucleotides), the slope of the repair time function is relatively large as compared to larger values of the ratio. Figure 3 captures the dynamics of only the repair part of the model, that is, it shows the repair times for a fixed amount of AP lesions.

Figure 3 

Stage 2 mathematical model of BER. The plot shows the time required to repair a fixed number of AP sites per 105 nucleotides with wild-type BER enzymes. The horizontal axis represents the starting number of AP sites per 105 nucleotides and the vertical axis the time for this number of sites to be repaired.

In the experiments addressed mathematically in Stage 1 of the model, we demonstrated that the CYP1B1-mediated oxidation of E₂ resulted in the formation of approximately 10 nmol/L 4-OHE₂-N7-Gua, which translates into 10 adducts per 10⁵ nucleotides per hour [2]. Since the formation of 4-OHE₂-N7-Gua leads to spontaneous depurination [13], we can assume an equal number of AP sites formed as a result of oxidative estrogen metabolism. This is one order of magnitude higher than the number of AP sites formed in normal leukocytes and fibroblasts, 10 AP sites per 10⁶ nucleotides per hour [13]. The living cell depends on efficient enzymatic repair of the continuous chemical damage inflicted on genomic DNA. Therefore, in the last stage of the model we combine the adduct formation model (Stage 1) and the adduct repair model (Stage 2). The mathematical model used in Stage 1 was chosen to fit batch data (a fixed amount of E₂). However, we consider the case when E₂ is continuously available and can vary over time and hence, there is continuous adduct formation and repair. In our modification of the Stage 1 model, we assume that the total amount of 2^′-deoxyguanosine is fixed, and it exists in two states: damaged (GUA_adduct) and undamaged or repaired (𝑦₉). The modified Stage 1 model, which accounts for the continuous nature of damage and repair, is then given by the following system of ordinary differential equations𝑑(4-OHE₂)=𝑘𝑑𝑡^cat1E_CYP1B1E₂𝐾_𝑚1+E₂−𝑉_maxQ4-OHE₂𝐾_𝑚Q+4-OHE₂,𝑑(E₂-3,4-Q)=𝑉𝑑𝑡^maxQ4-OHE₂𝐾_𝑚Q+4-OHE₂−𝑘_QE₂-3,4-𝑄,𝑑(GUA_adduct)𝑑𝑡=𝑘_𝑟(𝐸₂-3,4-Q)𝑦₇−𝑘_𝑠GUA_adduct. (3) The endpoint of the repair process is 𝑦₉(𝑡) which is the concentration of GUA(𝑡). The term, 𝑅(𝑡)=𝑘_𝑟𝑦₉(E₂-3,4-𝑄) where 𝑘_𝑟 is a constant, is the term that connects (2) and (3). In case the AP sites are not generated spontaneously but produced by DNA glycosylases, we can easily incorporate the enzymatic reaction catalyzed by glycosylases, which is not rate limiting for BER [9].

The combined Stage 3 model is composed of the differential in (2) and (3) along with their initial conditions. One can easily show that the system has a conservation lawGUA_adduct(𝑡)GUA₀1=1−GUA₀¹¹_𝑖=4𝑦_𝑖(𝑡). (4) Since this equation holds for all times, it calculates, at least in theory, the steady-state value of the unrepaired ratio.

3. Discussion

Numerous epidemiological studies have implicated estrogens in the development of breast cancer [1]. A pooled analysis of nine prospective studies of serum estrogen levels and breast cancer in over 2,400 postmenopausal women revealed a strong association of serum E₂ concentrations with breast cancer risk [14]. In Figure 4(a), we show the results of the combined model for a constant source of estrogen, E₂(𝑡)≡1, on the number of AP sites per 10⁵ nucleotides as well as the effect of doubling the E₂ level, E₂(𝑡)≡2. Clinical studies estimate that the doubling of serum E₂ levels confers a 1.3-fold increase in the risk of breast cancer [14]. Our model predicts a larger increase in the number of AP sites, indicating the involvement of other factors, such as enzyme activities. Enzymes involved in the production of DNA-damaging carcinogens and in repair of the resulting adducts are critical in maintaining the integrity of genomic DNA. Genetic variants of these enzymes, which occur in the general population, have been shown to play a role in altering specific reactions catalyzed by the individual enzymes. For example, we have shown that commonly occurring variants of CYP1B1 differ in their activity of producing 4-OHE₂ from the parent hormone E₂ [3, 15]. Similarly, the APE1 gene contains over 20 polymorphic sites, which have been assessed functionally as recombinant proteins or based on structural predictions [16, 17]. Several APE1 variants are associated with substantially reduced or increased activities. For example, L104R, E126D, and R237A exhibit 50% and D283G 90% reduction in repair capacity [16]. In view of these associations we used the model to examine the effect of variant enzyme haplotypes on the formation and repair of DNA adducts. In Figure 4(b), we illustrate the effect of replacing wild-type CYP1B1 with a variant having a 30% greater 𝑘_cat and wild-type APE1 with a variant having a 30% lower 𝑘_cat. As expected, the more active CYP1B1 variant increases DNA damage, and the less active APE1 variant results in decreased repair, which is reflected in higher levels of AP lesions for each. Their combined occurrence has an additive effect on AP lesions. Figure 4(c) shows the effect of different variants of CYP1B1 and APE1, which cause a decrease in number of AP sites relative to the wild-type enzymes. The panels in Figure 4 demonstrate the number the DNA adducts, which are constantly being created and repaired as modeled in (2) and (3), approaching steady values in the integrated process of adduct formation and repair.

(a)

(b)

(c)

Figure 4 

Effects of different concentrations of estrogen and variant haplotypes of CYP1B1 and APE1 on the formation and repair of DNA adducts. (a) Shows the effect of doubling the concentration of E2 on the number of AP lesions per 105 nucleotides over time using E2=1 (red curve) and E2=2 (blue curve). (b) Shows the effect of variant haplotypes of CYP1B1 and APE1 leading to an increase in the number of AP lesions per 105 nucleotides relative to wild-type enzymes. Each curve uses different assumptions about the kinetic constant, 𝑘cat, for the CYP1B1 and APE1 enzymes. The red curve depicts the wild-types with 𝑘cat1=1.6 (CYP1B1) and 𝑘cat5=192.168 (APE1). The blue curve uses a 30% increase in 𝑘cat1, keeping 𝑘cat5=192.168. The green curve employs a 30% reduction in 𝑘cat1, keeping 𝑘cat5=192.168. Finally, the purple curve uses both a 30% increase in 𝑘cat1 and a 30% reduction in 𝑘cat5. (c) Shows the effect of different variant haplotypes of CYP1B1 and APE1 that cause a decrease in the number of AP lesions per 105 nucleotides relative to wild-type enzymes. The wild-type is shown in red, a 30% decrease of CYP1B1 activity with APE1 activity held at its wild-type level in blue, 30% increase of APE1 activity with CYP1B1 activity held at its wild-type level in green, and a 30% decrease in CYP1B1 activity with a 30% increase of APE1 activity in purple.

A limitation of the model is its focus on BER and the omission of other DNA repair pathways. While the majority of adducts derived from 4-OHE₂ and E₂-3, 4-Q are depurinating N7 guanine adducts susceptible to BER, other stable, bulky DNA adducts that do not depurinate are also formed. However, the level of the stable adducts is three to four orders of magnitude lower than the level of depurinating adducts justifying the emphasis on BER [13, 18]. A vast number of depurination events occur under normal conditions involving not just guanine but other bases as well [19]. The importance of depurinating adducts derived from oxidative estrogen metabolism relative to the total spectrum of depurination events is presently unknown.

A critical strength of the model is that it can readily incorporate additional enzymes involved in the metabolism, as we have shown for CYP1A1, catechol-O-methyl transferase, and glutathione S-transferase P1 in the estrogen metabolism pathway [20]. Moreover, the model incorporates both phenotypic measures of estrogen exposure and genotypic data. This is schematically shown in Figure 5, which displays the interaction of estrogens, enzyme genotypes, and resulting AP sites per 10⁵ nucleotides as a three-dimensional graph. In designing Figure 5, we assumed that the difference in estrogen exposure between individual women is no more than twofold, with the tertiles 1.0, 1.5, and 2.0 used for the 𝑥-axis. This range is conservative since up to fivefold differences in breast tissue concentrations have been reported in [21]. Regardless of the scale used for the 𝑥-axis, the production of AP sites would be expected to be greater in women with more endogenous (more ovulatory cycles) or exogenous (hormone replacement therapy, oral contraceptives) estrogen exposure. The effect of wild-type and variant enzyme haplotypes on AP sites is displayed on the 𝑦-axis. It is evident from Figure 5 that the combined phenotypic and genotypic factors exert not just an additive but multiplicative effect on AP site levels, which are shown on the 𝑧-axis.

Figure 5 

Three-dimensional graph displaying the interaction of estrogen exposure, enzyme genotypes, and resulting AP sites per 105 nucleotides. We display estrogen exposure in tertiles on the 𝑥-axis. Estrogen exposure can be represented by actual E2 values, measured in pmol/L, in combination with semiquantitative estimates of each woman’s overall exposure to estrogen. The latter can be derived by taking into account her total years of ovulation as a function of current age, age at menarche, age at menopause, numbers of full-term pregnancies, and the dosage and duration of the use of exogenous estrogens. Thus, as far as the model is concerned, exogenous and endogenous estrogens can be combined although their precise contribution to estrogen exposure is presently unknown. The 𝑦-axis represents the combined effects of wild-type and variant enzyme haplotypes in the oxidative estrogen metabolism and BER pathways on AP levels. In theory, all enzyme genotype combinations shown in Figures 4(b) and 4(c) could be plotted. However for clarity, we have plotted only the wild-type and the lowest and highest variant haplotypes, separated into tertiles based on their respective AP sites per 105 nucleotides, which is represented on the 𝑧-axis. (The authors acknowledge Eric Parl for the design of this figure.)

In conclusion, we modeled the dynamic interaction between the estrogen-mediated DNA damage pathway and the DNA BER pathway in order to assess the overall impact on carcinogenesis. The model encompasses a set of differential equations representing the sequence of enzymatic reactions in both damage and repair pathways. The model allows us to examine the effect of phenotypic and genotypic factors such as different concentrations of estrogen and variant enzyme haplotypes on the formation and repair of DNA adducts. As a practical application, the model can be used to explore genetic factors in the damage-repair cycle (through the rate constants and estrogen levels) as a first step in understanding why some women develop breast cancer while others with similar circumstances do not. As better information about the kinetics in the model (e.g., rate constants for different haplotypes) becomes available, predictions about breast cancer risk will improve. Thus, the model may be useful for the construction of breast cancer risk models.

Acknowledgments

This work was supported by NIH (U54CA113007) and the Vanderbilt Integrative Cancer Biology Center.