A mathematical model for DNA damage and repair.

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 17beta-estradiol (E(2)) 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 E2 to 4-OHE2 and the estrogen quinone (E2-3,4-Q, Figure 1). The highly reactive quinone nonenzymatically attacks DNA and forms covalent adducts with bases, such as 4-OHE2-N7-Gua (GUAadduct), 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 E2 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, E2 denotes 17β-estradiol, 4-OHE2 the catechol estrogen, E2-3,4-Q the estrogen quinone, GUA the deoxyribonucleoside 2′-deoxyguanosine, and GUAadduct the estrogen adduct 4-OHE2-N7-Gua over time t. We assume that the starting points for the dynamic model are the initial concentration of GUA(0) = GUA0 and a prescribed time series is given for E2.

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 are

In the first reaction, the production of 4-OHE2 from E2 is an enzymatic reaction mediated by CYP1B1 that follows Michaelis-Menten kinetics with parameters K and kcat [3]. In the second reaction, E2-3,4-Q is produced from 4-OHE2 according to Hill-type kinetics. It has been shown that E2-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 kQ. 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 E2 leads to the formation of 4-OHE2 and deoxyribonucleoside adducts, such as 4-OHE2-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 E2, 4-OHE2, and GUAadduct.

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. S(t) denotes the concentration of DNA adducts. We define the source and removal function for the repair cycle functions, S(t) and R(t), as S(t) = kGUAadduct and R(t) = k(E2-3,4-Q)y9. Here, k is a modeling parameter expressed in [1/min ] which reflects the rate at which GUAadduct enters the repair cycle. The parameter k is the rate at which y9(t) interacts with E2-3, 4-Q to form GUAadduct. The two parameters, k and k, can be used to control the rates of entering and exiting of adducts and repaired lesions, respectively. For the simulations, we have chosen, k = k = 1,

There are several Michaelis-Menten parameters, kcat and K, for which experimental values have been estimated or measured (Table 1) [6]. The enzyme levels are denoted by EAPE1, EFEN1, EPol, EPol, ELig1, and ELig3. The variable, y9(t), represents the repaired GUAadduct, and if we set its initial condition in the system of differential equations as y9(0) = GUA0, it has the same meaning as GUA(t). 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 2000 nM [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].

Reaction subscript (i)	kcati [1/min]	Kmi [nmol/L]
5	192.168	34.7
6	4.500	500
7	49.020	210
8,11	1.278	56.7
9	36.000	100
10	8.130	39

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 (AP  lesions/105  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 E2 resulted in the formation of approximately 10 nmol/L 4-OHE2-N7-Gua, which translates into 10 adducts per 105 nucleotides per hour [2]. Since the formation of 4-OHE2-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 106 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 E2). However, we consider the case when E2 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 (GUAadduct) and undamaged or repaired (y9). 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

The endpoint of the repair process is y9(t) which is the concentration of GUA(t). The term, R(t) = ky9(E2-3,4-Q) where k 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 law

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 E2 concentrations with breast cancer risk [14]. In Figure 4(a), we show the results of the combined model for a constant source of estrogen, E2(t) ≡ 1, on the number of AP sites per 105 nucleotides as well as the effect of doubling the E2 level, E2(t) ≡ 2. Clinical studies estimate that the doubling of serum E2 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-OHE2 from the parent hormone E2 [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 kcat and wild-type APE1 with a variant having a 30% lower kcat. 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.

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, kcat, for the CYP1B1 and APE1 enzymes. The red curve depicts the wild-types with kcat = 1.6 (CYP1B1) and kcat = 192.168 (APE1). The blue curve uses a 30% increase in kcat, keeping kcat = 192.168. The green curve employs a 30% reduction in kcat, keeping kcat = 192.168. Finally, the purple curve uses both a 30% increase in kcat and a 30% reduction in kcat. (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-OHE2 and E2-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 105 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 x-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 x-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 y-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 z-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 x-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 y-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 z-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.