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 17β-estradiol (E2) 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.


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 2 to 4-OHE 2 and the estrogen quinone (E 2 -3,4-Q, Figure 1). The highly reactive quinone nonenzymatically attacks DNA and forms covalent adducts with bases, such as 4-OHE 2 -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.
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 estrogendeoxyribonucleoside 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 2 and 2 -deoxyguanosine, determine the extent of adduct formation and the correction by BER required to preserve the integrity of DNA.

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 2 denotes 17β-estradiol, 4-OHE 2 the catechol estrogen, E 2 -3,4-Q the estrogen quinone, GUA the deoxyribonucleoside 2 -deoxyguanosine, and GUA adduct the estrogen adduct 4-OHE 2 -N7-Gua over time t. We assume that the starting  In each graph the dots represent the experimental data [2], and the curves are derived from the mathematical model.
points for the dynamic model are the initial concentration of GUA(0) = GUA 0 and a prescribed time series is given for E 2 .
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-OHE 2 from E 2 is an enzymatic reaction mediated by CYP1B1 that follows Michaelis-Menten kinetics with parameters K m1 and k cat1 [3].
In the second reaction, E 2 -3,4-Q is produced from 4-OHE 2 according to Hill-type kinetics. It has been shown that E 2 -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 k 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 2 leads to the formation of 4-OHE 2 and deoxyribonucleoside adducts, such as 4-OHE 2 -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 2 , 4-OHE 2 , and GUA adduct . 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) = k s GUA adduct and R(t) = k r (E 2 -3, 4-Q)y 9 . Here, k s is a modeling parameter expressed in [1/ min] which reflects the rate at which GUA adduct enters the repair cycle. The parameter k r is the rate at which y 9 (t) interacts with E 2 -3, 4-Q to form GUA adduct . The two parameters, k s and k r , can be used to control the rates of entering and exiting of adducts and There are several Michaelis-Menten parameters, k cati and K mi , 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, y 9 (t), represents the repaired GUA adduct , and if we set its initial condition in the system of differential equations as y 9 (0) = GUA 0 , 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]. 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/10 5 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.
In the experiments addressed mathematically in Stage 1 of the model, we demonstrated that the CYP1B1mediated oxidation of E 2 resulted in the formation of approximately 10 nmol/L 4-OHE 2 -N7-Gua, which translates into 10 adducts per 10 5 nucleotides per hour [2]. Since the formation of 4-OHE 2 -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 6 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 2 ). However, we consider the case when E 2 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 (y 9 ). 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 Journal of Nucleic Acids Each curve uses different assumptions about the kinetic constant, k cat , for the CYP1B1 and APE1 enzymes. The red curve depicts the wild-types with k cat1 = 1.6 (CYP1B1) and k cat5 = 192.168 (APE1). The blue curve uses a 30% increase in k cat1 , keeping k cat5 = 192.168. The green curve employs a 30% reduction in k cat1 , keeping k cat5 = 192.168. Finally, the purple curve uses both a 30% increase in k cat1 and a 30% reduction in k cat5 . (c) Shows the effect of different variant haplotypes of CYP1B1 and APE1 that cause a decrease in the number of AP lesions per 10 5 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.
The endpoint of the repair process is y 9 (t) which is the concentration of GUA(t). The term, R(t) = k r y 9 (E 2 -3, 4-Q) where k r 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.

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 2 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 2 (t) ≡ 1, on the number of AP sites per 10 5 nucleotides as well as the effect of doubling the E 2 level, E 2 (t) ≡ 2.
Clinical studies estimate that the doubling of serum E 2 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 6 Journal of Nucleic Acids  Figure 5: Three-dimensional graph displaying the interaction of estrogen exposure, enzyme genotypes, and resulting AP sites per 10 5 nucleotides. We display estrogen exposure in tertiles on the x-axis. Estrogen exposure can be represented by actual E 2 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 wildtype 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 10 5 nucleotides, which is represented on the z-axis. (The authors acknowledge Eric Parl for the design of this figure.) 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 2 from the parent hormone E 2 [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 k cat and wild-type APE1 with a variant having a 30% lower k 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 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 2 and E 2 -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 5 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 xaxis. 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.
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.