In silico Modeling of Intracellular Diffusion and Reaction of Benzo[a]pyrene Diol Epoxide

(1)

In silico Modeling of Intracellular Diffusion and Reaction of Benzo[a]pyrene Diol Epoxide

Dreij K¹, Chaudhry QA², Jernström B¹, Hanke M² and Morgenstern R¹

1. Institute of Environmental Medicine, Karolinska Institutet, Stockholm, Sweden, 2. School of Computer Science and Communication, Royal Institute of Technology,

Stockholm, Sweden

March 23, 2012

Abstract

Several studies has suggested that glutathione conjugation of polycyclic aromatic hydrocarbons (PAHs) catalyzed by glutathione transferases (GSTs) are important factors in protecting cells against toxicity and DNA damage derived from these compounds. To further characterize the intracellular dynamics of PAH DEs and the role of GSTs in protection against DNA damage, we recently developed a PDE model using techniques for mathematical homogenization (Dreij K et al. PLoS One 6(8), 2011). In this study, we wanted to further develop our model by benchmarking against results from four V79 cell lines; control cells and cells overexpressing human GSTs A1- 1, M1-1 and P1-1. We used an approach of global optimization of the parameters describing the diffusion and reaction of the ultimate carcinogenic PAH metabolite benzo[a]pyrene diol epoxide to fit measured values from the four V79 cell lines. Numerical results concerning the formation of glutathione conjugates and hydrolysis were in good agreement with results from measurements in V79 cell culture. Cellular results showed significant protection by GST expression against formation of DNA adducts with more than 10- fold reduced levels compared to control cells. Results from the model using globally optimized parameters showed that the model cannot predict the protective effects of GSTs. Extending the model to also include effects from protein interactions and GST localization showed the same discrepancy. In summary, the results show that we have an incomplete understanding of the intracellular dynamics of the interaction between BPDE and GST that warrants further investigation and development of the model.

(2)

Introduction

When mammalian cells are exposed to foreign and potentially harmful compounds, a series of events takes place. Following uptake the substance is distributed in different intracellular compartments usually by passive diffusion, absorption and desorption. Depending on the lipophilicity, the majority of the compound is either dissolved in the aqueous phase, i.e. the cytoplasm, or in the lipophilic phase, i.e. the membranes. Simultaneously, biotransformation is catalyzed by different soluble and membrane bound enzymes with the purpose to render the substance non-toxic and suitable for excretion. The lipophilicity, which governs the partitioning between cytoplasm and membrane, is an important parameter in deter- mining bioavailability and thus toxic effects of a chemical. This has been shown in both biological and mathematical models for several groups of compounds including alkanols, chlorinated hydrocarbons and polycyclic aromatic hydrocarbons (PAHs) [1, 2].

PAHs are highly lipophilic widespread environmental pollutants and most probably carcinogenic to humans [3]. Probably the most well-investigated carcinogenic PAH is benzo[a]pyrene (BP). To exert its biological activity BP requires metabolic activation and subsequent DNA-adduct formation. One route of activation common to all studied PAHs includes the sequential action of cytochrome P450 and microsomal epoxide hydrolase and results in formation of bay or fjord region diol epoxides (DEs) [4]. The PAH DEs have been identified as the ultimate carcinogenic metabolites and DNA-adducts derived from DEs have shown to be mutagenic and carcinogenic [5].

In previous studies, to characterize the intracellular behavior of PAH DEs, we have examined the kinetics and role of soluble glutathione transferase (GST) catalyzed glutathione (GSH)-conjugation in protecting against DNA damage derived from several PAH DEs using in vitro [6–8] and mammalian cell models [9, 10].

Our data, together with others [11–13], implicate a central role for the GSTs Al- pha, Mu and Pi in protecting the genome against DE-induced damage and subsequent mutagenicity. To allow more sophisticated predictions of cellular behavior and thus toxicity of foreign compounds we recently developed a mathematical model describing the uptake and intracellular diffusion and reaction of PAH DEs [14]. This is the first model which includes the cytoplasmic membranes in such a diffusion reaction model, thus making it possible to study the effects of partitioning.

Here, we extend the framework of our previously developed model for further benchmarking against different GSTs. The results show that we have an incomplete understanding of the protective effects of GTSs which are much more pronounced than can be expected from their catalytic role and cytosolic location.

(3)

Material and Methods

Experimental Methods

Chemicals

Synthesis of the optically active (+)-anti-DE of benzo[a]pyrene has been performed according to literature methods. Additional chemicals and reagents were from Sigma-Aldrich and Fisher Scientific.

Cells and Treatment

The production and characterization of V79 Chinese hamster lung cells stably expressing human GSTs A1-1, M1-1 and P1-1 Val105, Ala114 and corresponding control cells has been previously described [10]. V79 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM), supplemented with 10% fe- tal bovine serum , 100 U/mL penicillin, 100mg/mL streptomycin, 1 mM sodium pyruvate (all from Gibco, Paisley, UK) and 0.1 mg/ml Hygromycin B (Invitrogen) at 37°C and 5% CO2 in a humidified environment. Cells were treated with DE as described [10]. In brief, cells in exposure medium (PBS containing 1 g/L glucose, pH 7.5) were incubated with 1 µM BPDE for up to 30 min. The extent of hydrolysis of BPDE during storage was determined to be in average 5%. At certain time points unreacted DE, GSH conjugates, tetrols and DNA adducts were analyzed by HPLC as previously described [6, 10]. Some samples were spiked with internal standard during the processing to adjust for loss of recovery.

Estimation of Protein Covalent Binding Constant

Confluent V79 cells were washed with PBS, pelleted at 1500 rpm for 5 min and washed twice with PBS. Cells were resuspended in 20 mM Tris-HCl buffer pH 7.5, 0.25 M sucrose, 0.5 mM EDTA, and sonicated 2 × 20 sec on ice. The cell lysate was spun 10 min at 4,000 rpm at 4°C and the subsequent supernatant 60 min at 40,000 rpm at 4°C. The resulting supernatant containing cytosolic proteins was concentrated by Amicon Ultra-4 devices (Millipore Corp, Billerica, MA, USA) 30 min at 3,750 rpm at 4°C. BPDE was added to cytosolic protein solution (0.05 to 5.0 mg/ml) to a final concentration of 20 µM and incubated at 37°C. At certain time points the reaction was stopped with the addition of alkaline 2-mercaptoethanol. The protein was precipitated with 4 × volume ice-cold acetone and the samples spun 10 min at 3,200 rpm at 4°C. The resulting pellet was washed three times with ice-cold acetone, dried under N2 (g) and subsequently dissolved in 5% aqueous SDS. Extent of BPDE-protein binding was analyzed by UV spectrophotometry, ε_345.7nm = 21000M⁻¹cm⁻¹.

(4)

Estimation of DNA Covalent Binding Constant

The reactivity of BPDE against DNA was determined using isolated nuclei from rat liver cells. Livers were removed from male Sprague Dawley rats (ca 180 g body weight) and the nuclear fraction isolated [15]. DNA was quantified using Hoechst 33258 staining as described [16]. Isolated nuclei corresponding to 100 µg DNA were incubated with 10 nmol BPDE in 100 µl 50 mM Tris-HCl buffer pH 7.5 at 37°C for up to 30 min in duplicates. At certain time points unreacted BPDE was trapped by addition of alkaline 2-mercaptoethanol and the DNA isolated by phenol chloroform extraction. DNA adducts were analyzed as tetrols released following hydrolysis in dilute HCl as described [17,18]. The tetrols were analyzed by HPLC using a Nova Pak 4 µm C18 (3.9 × 150 mm analytical column, Waters Inc.) and an elution system composed of 55% methanol in water. The tetrols were identified and quantitated by comparison with authentic standards obtained by acidic hydrolysis of pure DE enantiomers.

Model Implementation

Mathematical Model

The reactions considered in the model are shown in Figure 1. The model is ba- sically the same as that described previously [14] but has been extended to include the reaction of DE with cellular proteins. The system comprises three sub-domains; extracellular medium, cytoplasm and nucleus, separated by cellular and nuclear membranes, respectively as shown in Figure 1. The DEs, referred to as C, are taken up and equilibrates through the cellular sub-domains by diffusion. Because of the lipophilic nature of the modeled compound and its metabolites a major part of the molecules will be absorbed into the cellular membranes.

The partition coefficient, K_P, describes the equilibrium ratio of the concentration of C and its metabolites between any intracellular aqueous-lipid interface, i.e. membrane-cytosol. In the model, no reactions take place in the membranes.

The parameter values for the enzymatic and non-enzymatic reactions and diffusion coefficients are shown in Table 4. The geometric constants for the model, including volume of the cells and the medium and relative volume/thickness of the sub-domains/membranes were set as previously described [14].

Modeling Assumptions The following assumptions were earlier made in [14]

and discussed in detail there.

A1 On the microscopic view, we observe that the structure of the cytoplasm is very complex. It contains many thin membranes. We assume these membranes as layered structures.

(5)

A2 On the large scale, we assume that the volume of the cytoplasm contains an unordered set of small-scale substructures, which are uniformly distributed over the volume.

A3 The cell has spherical symmetry. With this assumption, the three dimensional model has been reduced to one dimensional computational model.

A4 The physical and chemical properties of the cytoplasm and the membranes are uniform. The diffusion coefficients will be constant in their respective sub-domains as the consequence of this assumption.

A5 We adopt the continuum hypothesis meaning that; the distribution of the substances will be described using the concentrations.

A6 The jump condition in the concentration between the two sub-domains can be conveniently described by the partition coefficients.

Using assumption A5, the distribution of the substances is described using concentrations. In order to distinguish among the concentrations within the different sub-domains, we add an index i.e. the quantities valid in sub-domain i carry the index i. In addition to the model given in [14], new reaction for protein binding has been added in cytoplasm and nucleus.

Detailed Equations The reactions and diffusion in Figure 1, give rise to the following system of reaction-diffusion partial differential equations.

• Sub-domain 1 (extracellular medium)

∂

∂ t

C₁= ∇ · (D1∇C₁) − k_UC₁,

∂

∂ t

U₁= ∇ · (D1∇U₁) + k_UC₁.

• Sub-domains 2 and 4 (cellular and nuclear membranes)

There is only diffusion process in this domain, and no chemical reaction is taking place.

For i = 2, 4, it holds:

∂

∂ tC_i= ∇ · (Di∇C_i),

∂

∂ tU_i= ∇ · (Di∇U_i).

(6)

• Sub-domain 3 (cytoplasm) σ_C,h

∂

∂ t

C₃= ∇ · (D3,C,h∇C₃) − (k_U,h+ k_B,h+ k_E,h)C₃, σ_U,h

∂

∂ tU₃= ∇ · (D3,U,h∇C3) + k_U,hC₃,

∂

∂ tB₃= k_B,hC₃,

∂

∂ tE₃= k_E,hC₃.

where σ_S,h,S = C,U represents the time scaling coefficient which is obtained due to the homogenization procedure, and is defined by,

σ_S,h =

(p_w+ p_l/K_P,S, x∈ acqeous part of cytoplasm 1/K_P_,S, x∈ lipid part of cytoplasm

where p_w and p_l denote the volume fraction of the aqueous and lipid parts, respectively, where it holds that p_w+ p_l= 1. K_P_,S represents the partition coefficient. ’h’ in the subscript of diffusion and rate constants stands for the homogenized diffusion and rate constants. The detailed derivation of the homogenized system, including the derivation of effective diffusion coefficients, D_3,S,h, and effective chemical constants in cytoplasm, can be found in [14].

• Sub-domain 5 (nucleus)

∂

∂ tC₅= ∇ · (D₅∇C₅) − (k_U+ k_A+ k_E)C₅,

∂

∂ tU₅= ∇ · (D5∇U₅) + k_UC₅,

∂

∂ tA₅= k_AC₅,

∂

∂ t

E₅= k_EC₅.

The basic values of all the chemical parameters have been summarized in Table 4.

Compartment Model Even after homogenization, the computational complex- ity for the model described above is rather complex. Therefore, we used a compartment model for some more expensive calculations. The complete set of equations for the compartment model is given below. In addition to to the model developed in [14], new reaction for protein binding has been added in the cytoplasm

(7)

and nucleus in the present model.

V₁d

dtC₁= DA₁ K_p,Cδ(C₃

σ_C−C₁) −V₁k_UC₁, V₁d

dtU₁= DA₁ K_p,Uδ(U₃

σ_U −U₁) +V₁k_UC₁, V₃d

dtC₃= DA₁

K_p,Cδ(C₁−C₃

σ_C) + DA₂

K_p,Cδ(C₅−C₃

σ_C) −V₃k_U,h+ k_B,h+ k_E,h σ_C C₃, V₃d

dtU₃= DA₁

K_p,Uδ(U₁−U₃ σU

) + DA₂

K_p,Uδ(U₅−U₃ σU

) +V₃k_U,h σ_C

C₃, V₃d

dtB₃= V₃k_B,h σ_C C₃, V₃d

dtE₃= V₃k_E,h σ_C

C₃, V₅d

dtC₅= DA₂ K_p,Cδ(C₃

σC

−C₅) −V₅(k_U+ k_A+ k_E)C₅, V₅d

dtU₅= DA₂ K_p,Uδ(U₃

σ_U −U₅) +V₅k_UC₅, V₅d

dtA₅= V₅k_AC₅, V₅d

dtE₅= V₅k_EC₅.

Dimensional Analysis For the purpose of dimensional analysis, a typical length (X ), time (τ), and concentration ( ˜C) have been defined. We choose them as follows:

• Length scale: Radius of a cell. It is computed from the volume of a cell under the assumption that the cell has a spherical shape. If we denote the radius of cell by R, then X = R.

• Time scale: We take the rate constant of U in the extracellular medium as gauge value. Hence, τ = k_U.

• Concentration: The choice is not critical since the system is linear in all concentrations. We choose ˜C = C₀, the initial concentration of C in the extracellular medium.

(8)

Using the new scale variables, the compartment equations take the form, d

dτC˜₁= D_1,C,V₁(C˜₃

σ_C− ˜C₁) − ˜C₁ d

dτU˜₁= D_1,U,V₁(U˜₃

σ_U − ˜U₁) + ˜C₁ d

dτC˜₃= D_1,C,V₃( ˜C₁−C˜₃

σ_C) + D_2,C,V₃( ˜C₅−C˜₃

σ_C) − (K_U+ K_B+ K_E) ˜C₃ d

dτU˜₃= D_1,U,V₃( ˜U₁−U˜₃

σ_C) + D_2,U,V₃( ˜U₅−U˜₃

σ_U) + K_UC˜₃ d

dτB˜₃= K_BC˜₃ d

dτE˜₃= K_EC˜₃ d

dτC˜₅= D_2,C,V₅(C˜₃

σ_C− ˜C₅) − (˜k_U+ ˜k_A+ ˜k_E) ˜C₅ d

dτU˜₅= D_2,U,V₅(U˜₃ σU

− ˜U₅) + ˜k_UC˜₅ d

dτA˜₅= ˜k_AC˜₅ d

dτ

E˜₅= ˜k_EC˜₅

where for S = C,U, D_1,S,V₁ = _˜^{D ˜}^˜^A¹

V1Kp,Sδ˜, D_i,S,V₃ = _˜^{D ˜}^˜^Aⁱ

V3Kp,Sδ˜ (for i = 1, 2), D_2,S,V₅ =

D ˜˜A₂

V˜₅K_p,Sδ˜, K_U = ^˜k^U,h

σC , K_B= ^˜k^B,h

σC , K_E = ^˜k^E,h

σC .

Table 1 and 2 provide an overview of the scaled parameters. From Table 2, we see that the diffusion related parameters have very small scaled values. Thus, their time scale is beyond the resolution of the measurements.

Global Parameter Fitting

The parameters have been optimized to fit the measured values. The quality of a fit is measured by the total error between the numerical results of the model and the experimental data. The usual way of minimizing this error is done by defining the objective function,

f(k, y) =1 2

n

∑

i=1

w_i(ϕi(k) − y_i)²

(9)

where ϕ represents the vector of numerical results depending upon the set of parameters (k), and y represents the experimental observations. wi stands for the weights.

Sensitivity measures the changes in the estimated parameters (here, k) against perturbations in the data y. It is given by,

s(k, y) = dk dy

where k = k(y) is given by the solution of the minimization problem min

k

f(k, y).

The weights w_iare determined as follows:

1. The measured and computed values have been normalized using represen- tative sizes for their concentrations.

2. The reliability of the measured values has been taken into account by mod- ifying the weights by the scaled standard deviations of the individual measurements.

3. In a first series of optimization, the weights have been used as described above giving roughly equal influence to all measurements. In a second numerical experiment, the weight for A has been increased whereas the weights for other species were decreased. The aim of doing this is to determine optimized parameters to accurately predict the protection of GST against DNA adducts in all 4 cell lines.

This provides us with

W = S_d⁻¹S⁻¹

where S_d= diag( ¯s₁, ¯s₂, ..., ¯s_n) is the diagonal matrix of scaling values, while S = diag(e₁, e₂, ..., e_n) contains the standard deviations of the measured values.

Finally, we obtain

s(k, y) = dk

dy = H⁻¹ϕ⁰(k)^TW, (1)

where H denotes the Hessian of ϕ(k).

For a qualitative estimation of the influence of different parameters on the mathematical model, the compartment model was used. The PDE model could not be used for the purposes of this sensitivity analysis because it is computationally too complex and the necessary gradient information is not available.

(10)

The parameter fitting problem is formulated in terms of a constrained optimization problem where the constraints consist essentially of non-negativity con- ditions. We observed that, within the search region, there are many local minima. Therefore, optimization calls for global technique such as theMultistart method [19] from Matlab’s Optimization Toolbox, which has the ability to per- form a global search. A similar kind of approach has been followed in [20].

The Multistart method uses gradient-based methods to return local and global minima. It runs a local solver from all starting points. It also provides the users the flexibility in choosing different local nonlinear solvers. For our problem, we have usedfmincon as local solver.

fmincon is a gradient based optimization function available in Matlab’s Op- timization Toolbox [21]. It minimizes the function with several variables, taking into account their available bounds as well. In order to use the trust-region- reflective algorithm, one needs to provide the gradient function. For the compartment model, it consists essentially of the solution of the equation of first variation for the system of differential equations. Thus, it is easily available.

For the estimation of the Hessian matrix, we numerically approximated the Jacobian of the gradient function using numjac provided by Matlab [22]. Tech- nically it was not possible to use theMultistart method with the trust-region- reflective algorithm. Therefore,fmincon was used (without Multistart) with the trust-region-reflective algorithm for the sensitivity analysis, and with the initial guess obtained from theMultistart method.

For our optimization problem using compartment model, we included K_p,C, K_p,U, kU, k_A, kB_P, kB_A, kB_M and kB_C(parameters were included in the optimization code in that order). The vector of eigenvalues (λall) of the Hessian matrix is determined at the fitted parameters, resulting in

λall =







8.2841 × 10¹² 7.3807 × 10¹² 3.6507 × 10¹² 2.0616 × 10¹² 1.5728 × 10¹¹ 1.0514 × 10¹¹ 2.6299 × 10⁸ 1.9819 × 10⁻¹







whereas the matrix Vall of eigenvectors of the Hessian matrix is:

(11)

Vall =







−0.3046 0.0422 −0.0339 −0.6168 −0.2079 −0.1983 −0.6643 −0.0000 0.0000 0.0000 −0.0000 −0.0000 −0.0000 −0.0000 −0.0000 1.0000

−0.1037 0.0121 −0.0120 −0.1714 −0.7607 −0.3042 0.5370 0.0000

−0.0011 −0.0006 0.0006 0.0018 0.4430 −0.8877 0.1251 0.0000 0.0471 −0.0074 0.0149 0.7232 −0.4056 −0.2696 −0.4870 −0.0000 0.3695 −0.9135 −0.0107 −0.1333 −0.0665 −0.0431 −0.0695 −0.0000 0.8704 0.4044 −0.0173 −0.2185 −0.1126 −0.0738 −0.1124 −0.0000 0.0067 −0.0011 0.9990 −0.0391 −0.0131 −0.0076 −0.0116 −0.0000







From the vector of eigenvalues of Hessian matrix, we see that one eigenvalues is (almost) zero. Due to the zero eigenvalue, the fitted parameters cannot be uniquely determined. From the matrix of corresponding eigenvectors we expect that, we cannot determine Kp,U given the provided data. Thus, we fixed the parameter K_p,Uat its measured value which is given in the Table 4, whereas the other parameters were included in the optimization code for data fitting. Once again the eigenvalues of the Hessian matrix of the reduced system were obtained. The sensitivity of the optimal parameters was calculated using Eq. 1. The relative error in the fitted parameters are given in the Table 3.

The PDE model was implemented in COMSOL Multiphysics3.5a [23] and interfaced to Matlab’s Optimization Toolbox. The globally optimized parameters are found in Table 5.

Results and Discussion

Previous work from our laboratory and by others has shown that GST and GSH are important factors in protecting cells against toxicity and DNA damage derived from PAHs [6–13]. To further characterize the intracellular dynamics of PAH DEs and the role of GSTs in protection against DNA damage we recently developed a PDE model using techniques for mathematical homogenization [14]. Initial numerical results concerning reaction and metabolism of PAH DEs were in good agreement with results from cellular experiments. Furthermore, compared to a compartment model with well-stirred compartments our model allowed for a more realistic representation. In this study we wanted to further develop our model by benchmarking against results from four V79 cell lines; control cells and cells overexpressing human GSTs A1-1, M1-1 and P1-1. In order to investigate how the model mimics cellular exposure, uptake, metabolism and reaction of BPDE, data from in vitro experiments and cells in culture describing the partitioning, intracellular metabolism, and reactivity of BPDE were collected. The reaction and partition constants governing the dynamics can be found in Table 4. When measurements of the different species were performed the total amounts of BPDE, tetrols and GSH conjugates were analyzed. Accordingly, the model was utilized to mimic this situation by adding the intra and extra cellular profiles.

(12)

Modeling Glutathione Conjugation and Tetrol Formation

As can be seen in Figure 2, the model demonstrated reasonable agreement with the results of cellular measurements for the conversion of BPDE and formation of tetrols for all three V79 cell lines overexpressing GSTs as well as the control V79 cell line. However, a significant difference between model and cells in the formation of GSH conjugates for the GSTA1-1, GSTM1-1 and control cells was observed. While formation of GSH conjugates was well predicted for GSTP1-1 cells, the same reaction was overestimated by a factor of 2-3 for the two other GST overexpressing cell lines and underestimated by factor in the same range for the control cells. The discrepancy between model and control cells is explained by using an arbitrary low k_cat/K_M for the constitutive GST based on our previous studies showing that constitutive V79 GST essentially lacks activity against BPDE [24]. The overestimation of formation of GSH conjugates for GSTs A1-1 and M1-1 confirm that there is no exact correlation between activities obtained from experiments with purified enzymes and the activities observed in complex cellular systems [10].

In our earlier paper we proposed that the observed differences between PDE model and cells in part might be explained by the reaction of BPDE with other cellular macromolecules, such as proteins. To further investigate this the rate of binding of BPDE to cytosolic protein (k_E) was determined using cytosolic fractions from control V79 cells as described in material methods and estimated to be 4.3 × 10⁻⁴s⁻¹. In contrast to earlier studies performed in different cell- like systems showing that up to about 10% of the total reaction of BPDE could be accounted for as covalent binding to proteins the addition of our experimen- tally obtained protein binding constant had minimal impact on the modeled cellular diffusion and reaction of BPDE (resulting in up to 0.12% protein bound BPDE) [25, 26].

Next we used a global optimization technique,

Multistart

[19] to optimize the parameters included in the model that are accessible to optimization for the given set of measurements, to fit the cellular results. With this technique, the solver searches for the global optimum solution of the parameters. The resulting fitted parameters can be seen in Table 5 with the corresponding curves plotted in Figure 2 (dashed line). The estimation of the relative errors has been computed using the compartment model. As can be seen, most parameters were changed within one order or magnitude with the exception of rate of GSH conjugation in control cells which can be explained as discussed above.

(13)

Modeling DNA Adduct Formation

In order to include the formation of DNA adducts in the model to study the impact of GST metabolism, the rate of DNA binding of BPDE was determined using isolated nuclei from rat liver as described in material and methods. Using intact nuclei allowed for estimation of a more relevant reaction rate in relation to intracellular diffusion and reaction of BPDE compared to using naked DNA. The DNA binding reaction constant was estimated to 6.2 × 10⁻³s⁻¹. This value is 10-fold lower than corresponding binding rate to naked DNA [27], most likely reflecting the impact of chromatin structure and nuclear proteins on DNA binding. This difference is in agreement with studies on the distribution of adducts derived from BPDE in the eukaryotic genome demonstrating up to 5-fold higher frequency of adducts in linker and genetically active DNA than in nucleosomal and inactive DNA [28, 29].

Comparing levels of DNA adducts in the model with results from the 4 cell lines (taken from [10]) showed that the model could not predict the protection of GSTs against the formation of DNA adducts (Figure 3A). While cell results showed up to more than 10-fold reduction of DNA adducts in the presence of GST, the model only showed up to 15% reduction.

Results from the global optimization to fit the parameters to the results from all 4 cell lines showed the reverse discrepancy between model and cells. Although the model now could accurately predict level of DNA adducts in GST overexpressing cells, the level in control cells was underestimated. However, there was only data for A from one time point available, reducing the impact of this parameter (k_A) on the fitted values. Consequently, to be able to determine optimized parameters to accurately predict the protection of GST against DNA adducts in all 4 cell lines the fitting of k_A was given full weight (Figure 3A). This was done analytically by increasing the weights for A, whereas the weights for other species were decreased. The resulting parameters can be found in Table 5. While lipophilicity and reactivity against water and DNA were within one order of magnitude from both the original and globally fitted parameters, the rates of GSH conjugation (k_B) showed larger deviations from the starting parameters. According to this analysis, to accurately predict the protection of GSTs against DNA adducts the difference in rate of GSH conjugation between overexpressing and control cells needed to be larger. These results suggest that the model needs to be modified or extended in order to properly model the protection from DNA adduct formation in cells expressing GSTs.

(14)

Impact of Intracellular Interactions of BPDE on DNA Adduct Formation

Recent papers have shown that GSTs can have different localizations and suggested the importance of GST localization for its function as a phase II enzyme protecting against DNA damage. Stella et al showed that up to 10% of the cytosolic pool of GST Alpha class enzymes were localized in the nucleus in rat hepatocytes [30]. Furthermore, Fabrini et al recently showed that detoxifying and antioxidant enzymes, including GSTs, can form a nuclear shield resulting in an increased protection against chromosomal oxidative damage (26). Consequently, nuclear localization or a shield of GSTs could have an impact on the protection against DNA damage derived from BPDE. We tested this hypothesis by introducing these GST localizations into the PDE model. Data describing the width of the shield and the GST concentration was taken from [31]. Comparing control and GSTP overexpressing cells the results showed that placing up to 100% of the cytosolic GST protein pool in the nucleus does not affect the formation of DNA adducts, neither did the inclusion of a nuclear shield (Figure 3B).

Since BPDE is highly lipophilic a large proportion of the intracellular molecules will partition and diffuse inside the membranes throughout the cell. Consequently, the interaction of BPDE with membranes could affect the extent of GSH conjugation and reaction with for example DNA. To test if the structural organization of the membranes could affect the intracellular dynamics of BPDE and thus protection of GSTs against DNA adducts the orientation of the membranes were changed. In the homogenized PDE model the orientation of the membranes is also homogenized here, the orientation of the membranes was changed to be either fully parallel or transversal to the gradient of BPDE. This could easily be implemented by a slight modification of the homogenization algorithm. As can be seen in Figure 3B neither of these changes had an effect on DNA adduct levels.

The cytosol contains in the range of 300 mg/ml protein and earlier studies have suggested a role of intracellular PAH-protein interactions in the regulation of metabolism and could thus potentially also influence the extent of DNA adducts formed [32, 33]. To test this we introduced two different scenarios of BPDE- protein interactions in the model. The first representing a more unspecific interaction with cytoplasmic proteins using a dissociation constant (k_d) in the micromo- lar range (kon= 1 × 10⁹M⁻¹s⁻¹, k_off= 1 × 10³s⁻¹). As can be seen in Figure 3B, although this decreased the levels of A it could still not predict the GST protection against DNA damages. Previous studies regarding the intracellular dynamics of BP reported the presence of a BP binding protein with a k_d in the nanomolar range [34]. It was suggested that the GSTs could be responsible for this. This was tested by introducing a second more specific protein interaction with the GSTs (kon= 1 × 10⁹M⁻¹s⁻¹, k_off= 1s⁻¹) and comparing the results from modeling for-

(15)

mation of A in control and GSTP overexpressing cells. Although this lead to a small reduction in levels of A in the GSTP overexpressing cells compared to control cells it did not fit with the cellular results. Moreover, the addition of protein interactions significantly affected the kinetics of the rest of the metabolism evident as prolonged intracellular half-life of BPDE resulting in significantly decreased levels of intracellular tetrols in parallel to increased levels of GSH-conjugates no longer matching the cellular results (results not shown).

In summary, none of the different possible intracellular interactions of BPDE that was modeled here could explain the protective effects of GSTs. This shows that we have an incomplete understanding of the intracellular dynamics of the interaction between BPDE and GST and warrants further investigation. Interest- ingly, previous studies have suggested that the GSTs interact with cellular membranes possibly resulting in increased accessibility of substrates present in the membranes [30, 35]. This possibility will be tested in future simulations, which will however require extensive development of the model.

(16)

Figure 1.Schematic diagram showing the reactions and diffusion of BPDE in and around one cell.

(17)

Figure 2. Comparison between results from V79 cells, PDE model and globally optimized parameters.

Simulated amounts of the different species were generated by running the model for 30 min.

Parameters used for the PDE model and from global fit are found in Table 1 and 2, respectively. Results from cellular experiments show mean ± SEM, n=2.

(18)

Figure 3.Modeling the protection of GST against DNA damages.A, comparison between results from V79 cells, PDE model, global optimization and full weight for kA. B, Effects of changed parameters related to the intracellular dynamics and interactions of BPDE on the protective effects of GST.

(19)

Table 1: Problem parameters and their scaled values

Parameter Value Scaled

cell radius 8.947×10⁻⁶m 1

k_U 4.3×10⁻³s⁻¹ 1

C₀ 10⁻⁶M 1

k_E 4.3×10⁻⁴s⁻¹ 10⁻¹ k_B_P_,h^b 2.4291×10⁻¹s⁻¹ 5.649×10¹ k_B_A_,h^b 8.1168×10⁻²s⁻¹ 1.8876×10¹ k_B_M_,h^b 1.5667×10⁻¹s⁻¹ 3.6434×10¹ k_B_C_,h^b 5.9683×10⁻⁵s⁻¹ 1.3880×10⁻² k_A 6.2×10⁻³s⁻¹ 1.4419 k_E_,h 3.2079×10⁻⁴s⁻¹ 7.4603×10⁻²

D 10⁻¹²m²s⁻¹ 2.9052

k_U,h 1.0444×10⁻³s⁻¹ 7.4603×10⁻¹

K_p,C 1.2×10⁻³ 1.2×10⁻³

K_p,U 8.3×10⁻³ 8.3×10⁻³

σ_C 212.4 212.4

σ_U 31.3 31.3

bCalculated for GSTP cell line.

(20)

Table 2: Scaled values of diffusion and reaction related parameters.

Diffusion related parameters

Reaction related parameters

Measured data available at time scales

parameter scaled value parameter scaled value Species scaled value

D_1,C,V₁ 3.8376e-005 k_U 1 C^c 0,0.26,0.52,1.3,2.6,5.2,7.7

D_1,U,V₁ 2.6543e-004 K_U 284.69 U^c 0,0.26,0.52,1.3,2.6,5.2,7.7 D_1,C,V₃ 1.2919e-007 K_B_P 3.76 B₃_P 0,0.26,0.52,1.3,2.6,5.2,7.7 D_1,U,V₃ 8.9355e-007 K_B_A 11.25 B₃_A 0,1.3,2.6,5.2,7.7

D_2,C,V₃ 3.2423e-007 K_B_M 5.83 B₃_M 0,1.3,2.6,5.2,7.7

D_2,U,V₃ 2.2426e-006 K_B_C 15300 B₃_C 0,5.2

D_2,C,V₅ 1.0884e-007 k_A 0.69 A₅ 0,5.2

D_2,U,V₅ 7.5284e-007 k_E 10 E Not available

σ_C 212.4 K_E 2847

σ_U 31.3

cOnly sum in all sub-domains available, not available in indivisual sub-dmains.

(21)

Table 3: Fitted parameter obtained using PDE model and relative errors of the optimal parameters using CM.

Parameter Starting value Fitted value ^e(Rel.Er) (%) K_p,C 1.2×10⁻³ 2.1680×10⁻³ 0.46

K_p,U 8.3×10⁻³ Not fitted

k_U 4.3×10⁻³ 5.5708×10⁻³ 1.75 k_A 6.2×10⁻³ 1.2855×10⁻³ 18.09

dK_B_P 3.256×10⁻¹ 3.5923×10⁻¹ 5.08

dK_B_A 1.088×10⁻¹ 5.0950×10⁻² 6.06

dK_B_M 2.100×10⁻¹ 8.2393×10⁻² 4.16

dK_B_C 8×10⁻⁵ 8.5848×10⁻³ 17.64 k_E 4.3×10⁻⁴ Not fitted

D₂, D₄ 1×10⁻¹² Not fitted

dG∗ K_C

eRelative error (Rel.Er) was computed using CM.

(22)

Table 4: Constants describing diffusion and reaction of BPDE for the model.

Symbol Constant Value Ref.

D1 Diffusion coefficient in extracellular medium [m²s⁻¹] 10⁻⁹ [36, 37]

D₂, D₄ Diffusion coefficient in cell/nuclear membrane [m²s⁻¹] 10⁻¹² [38]

D_3,lt Diffusion coefficient in cytoplasm membranes/tangential [m²s⁻¹] 10⁻¹² [38]

D_3,ln Diffusion coefficient in cytoplasm membranes/normal [m²s⁻¹] 10⁻¹² [38]

D_3,w Diffusion coefficient in cytoplasm [m²s⁻¹] 2.5 × 10⁻¹⁰ a D₅ Diffusion coefficient in nucleus [m²s⁻¹] 2.5 × 10⁻¹⁰ a

c0 Initial concentration in extracellular medium [M] 10⁻⁶

K_P,C DE partition coefficient 1.2 × 10⁻³ b

K_P_,U Tetrol partition coefficient 8.3 × 10⁻³ b

k_U Solvolytic reactivity [s⁻¹] 4.3 × 10⁻³ [10]

k_A DNA reactivity [s⁻¹] 6.2 × 10⁻³ c

k_E Protein reactivity [s⁻¹] 4.3 × 10⁻⁴ c

kcat/K_M Catalytic efficiency [M⁻¹s⁻¹]

GSTA1-1 3.2 × 10³ [6]

GSTM1-1 10.5 × 10³ [7]

GSTP1-1 3.7 × 10³ [7]

Constitutive GST 25 d

G₁ Concentration of cellular GST – Phase II [M]

GSTA1-1 3.4 × 10⁻⁵ [10]

GSTM1-1 2.0 × 10⁻⁵ [10]

GSTP1-1 8.8 × 10⁻⁵ [10]

Constitutive GST 3.2 × 10⁻⁶ [10]

a Based on DH₂O≈ 10⁻⁹m²s⁻¹ of benzo[a]pyrene [36, 37] and the relationship that D_intracell≈ 2.5D_H₂_O[39, 40].

b Determined using ALOGPS 2.1 software [41].

c This study.

d arbitrary low number [24]

(23)

Table 5: Comparison of constants used in the PDE model, before and after global parameter optimization.

Parameter PDE Globally optimized Full weight k_A

K_P_,C 0.0012 0.0022 0.000737

k_U [s⁻¹] 0.0043 0.0056 0.005514

k_A[s⁻¹] 0.0062 0.0013 0.006679

k_B,GSTP[s⁻¹] 0.33 0.36 7.508407

k_B,GSTA[s⁻¹] 0.11 0.051 3.923624

k_B,GSTM[s⁻¹] 0.21 0.082 25.02053

k_B,Control[s⁻¹] 0.00008 0.0086 0.000463 k_B= k_cat/K_M× G, (see Table 1).

(24)

References

[1] S.B. Dorn, G.H. Degen, H.M. Bolt, J. van der Louw, F.A. van Acker, D.J.

van den Dobbelsteen, and J.P. Lommerse. Some molecular descriptors for non-specific chromosomal genotoxicity based on hydrophobic interactions.

Arch Toxicol., 82:333–338, 2008.

[2] J. Jacob. The Significance of Polycyclic Aromatic Hydrocarbons as Environ- mental Carcinogens. 35 Years Research on PAH - a Retrospective. Polycycl.

Aromat. Comp., 28:242–272, 2008.

[3] World Health Organization and IARC. Some non-heterocyclic polycyclic aromatic hydrocarbons and some related occupational exposures. IARC Press, Lyon, France, 2010.

[4] D.R. Thakker, H. Yagi, W. Levin, A.W. Wood, A.H. Conney, and D.M. Je- rina. Polycyclic aromatic hydrocarbons: Metabolic activation to ultimate carcinogens. In M.W. Anders, editor, Bioactivation of Foreign Compounds, pages 177–242. Academic Press, London, UK, 1985.

[5] R.G. Harvey. Polycyclic aromatic hydrocarbons: Chemistry and carcino- genicity. Cambridge University Press, Cambridge, UK, 1991.

[6] B. Jernström, M. Funk, H. Frank, B. Mannervik, and A. Seidel. Glutathione S-transferase A1-1-catalysed conjugation of bay and fjord region diol epoxides or polycyclic aromatic hydrocarbons with glutathione. Carcinogenesis, 17(7):1491–1498, 1996.

[7] K. Sundberg, M. Widersten, A. Seidel, B. Mannervik, and B. Jernström.

Glutathione conjugation of bay- and fjord-region diol epoxides of polycyclic aromatic hydrocarbons by glutathione transferases M1-1 and P1-1. Chem.

Res. Toxicol., 10(11):1221–1227, 1997.

[8] K. Dreij, K. Sundberg, A.S. Johansson, E. Nordling, A. Seidel, B. Persson, B. Mannervik, and B. Jernström. Catalytic activities of human alpha class glutathione transferases toward carcinogenic dibenzo[a,l]pyrene diol epoxides. Chem. Res. Toxicol., 15(6):825–831, 2002.

[9] K. Sundberg, A. Townsend, A. Seidel, and B. Jernström. Glutathione conjugation and DNA adduct formation of diol epoxides in V79 cells expressing human glutathione transferase P1-1. Polycycl. Arom. Comp., 21(1):123–133, 2000.

(25)

[10] K. Sundberg, K. Dreij, A. Seidel, and B. Jernström. Glutathione conjugation and DNA adduct formation of dibenzo[a,l]pyrene and benzo[a]pyrene diol epoxides in V79 cells stably expressing different human glutathione transferases. Chem. Res. Toxicol., 15(2):170–179, 2002.

[11] M.E. Kushman, S.L. Kabler, S. Ahmad, J. Doehmer, C.S. Morrow, and A.J. Townsend. Protective efficacy of hGSTM1-1 against B[a]P and (+)- or (-)-B[a]P-7,8-dihydrodiol cytotoxicity, mutagenicity, and macromolecular adducts in V79 cells coexpressing hCYP1A1. Toxicol. Sci., 99:51–57, 2007.

[12] M.E. Kushman, S.L. Kabler, M.H. Fleming, S. Ravoori, R. C. Gupta, J. Doehmer, C.S. Morrow, and A.J. Townsend. Expression of human glutathione S-transferase P1 confers resistance to benzo[a]pyrene or benzo[a]pyrene-7,8-dihydrodiol mutagenesis, macromolecular alkylation and formation of stable N2-Gua-BPDE adducts in stably transfected V79MZ cells co-expressing hCYP1A1. Carcinogenesis, 28:207–214, 2007.

[13] M.E. Kushman, S.L. Kabler, S. Ahmad, J. Doehmer, C. S. Morrow, and A.J.

Townsend. Cytotoxicity and mutagenicity of dibenzo[a,l]pyrene and (+/-)- dibenzo[a,l]pyrene-11,12-dihydrodiol in V79MZ cells co-expressing either hCYP1A1 or hCYP1B1 together with human glutathione-S-transferase A1.

Mutat. Res., 624:80–87, 2007.

[14] K. Dreij, Q. A. Chaudhry, B. Jernstrom, R. Morgenstern, and M. Hanke.

A method for efficient calculation of diffusion and reactions of lipophilic compounds in complex cell geometry. PLoS One, 6:e23128, 2011.

[15] M. Kurokawa, Y. Hibino, and N Sugano. Associations of benzo[a]pyrene with rat-liver chromatin and the chromatin protein. Chem. Biol. Interact., 39:17–30, 1982.

[16] T. Araki, A. Yamamoto, and M Yamada. Accurate determination of DNA content in single cell nuclei stained with Hoechst 33258 fluorochrome at high salt concentration. Histochemistry, 87:331–338, 1987.

[17] R.O. Rahn, S.S. Chang, J.M. Holland, and L.R. Shugart. A fluorometric- HPLC assay for quantitating the binding of benzo[a]pyrene metabolites to DNA. Biochem. Biophys. Res. Commun., 109:262–268, 1982.

[18] K. Alexandrov, M. Rojas, O. Geneste, M. Castegnaro, A. M. Camus, S. Petruzzelli, C. Giuntini, and H. Bartsch. An improved fluorometric assay for dosimetry of benzo(a)pyrene diol-epoxide-DNA adducts in smokers’

(26)

lung: comparisons with total bulky adducts and aryl hydrocarbon hydroxy- lase activity. Cancer Res., 52:6248–6253, 1992.

[19] The MathWorks Inc. Global Optimization Toolbox 3.1.1, 2011.

[20] Novac M, Vladu E, Novac O, and Gordan M. Aspects Regarding the Op- timization of the Induction Heating Process using Fmincon, Minimax func- tions and Simple Genetic Algorithm. JEEE, 2(2):64–69, 2009.

[21] The MathWorks Inc. Matlab Optimization Toolbox 6.0, 2011.

[22] The MathWorks Inc. MATLAB 7.12.0.635 (R2011a), 2011.

[23] COMSOL AB, Stockholm. Comsol Multiphysics 3.5a, 2010.

[24] K. Sundberg, B. Jernström, and S. Swedmark. Studies on the differential in- hibition of glutathione conjugate formation of (+)-anti-benzo[a]pyrene 7,8- dihydrodiol 9,10-epoxide and 1-chloro-2,4-dinitrobenzene in V79 Chinese hamster cells. Biochem. J., 349:693–696, 2000.

[25] N.B. Islam, D.L. Whalen, H. Yagi, and D.M. Jerina. Kinetic studies of the reactions of benzo[a]pyrene-7,8-diol 9,10-epoxides in aqueous solutions of human serum albumin and nonionic micelles. Chem. Res. Toxicol., 1:398–

402, 1988.

[26] B. Jernström, L. Dock, and M. Martinez. Metabolic activation of benzo[a]pyrene-7,8-dihydrodiol and benzo[a]pyrene-7,8-dihydrodiol- 9,10epoxide to protein-binding products and the inhibitory effect of glutathione and cysteine. Carcinogenesis, 5:199–204, 1984.

[27] N.E Geacintov, H. Yoshida, V. Ibanez, S.A. Jacobs, and R.G. Harvey. Con- formations of adducts and kinetics of binding to DNA of the optically pure enantiomers of anti-benzo(a)pyrene diol epoxide. Biochem. Biophys. Res.

Commun., 122:33–3–9, 1984.

[28] A. Kootstra, T.J. Slaga, and D.E. Olins. Interaction of benzo[alpha]pyrene diol-epoxide with nuclei and isolated chromatin. Chem. Biol. Interact., 28:225–236, 1979.

[29] F.O. Obi, A.J. Ryan, and M.A. Billett. Preferential binding of the carcinogen benzo[a]pyrene to DNA in active chromatin and the nuclear matrix. Carcino- genesis, 7:907–913, 1986.

(27)

[30] L Stella, V. Pallottini, S. Moreno, S. Leoni, F. De Maria, P. Turella, G. Fed- erici, R. Fabrini, K.F. Dawood, M.L. Bello, J.Z. Pedersen, and G. Ricci.

Electrostatic association of glutathione transferase to the nuclear membrane.

Evidence of an enzyme defense barrier at the nuclear envelope. J. Biol.

Chem., 282:6372–6379, 2007.

[31] R. Fabrini, A. Bocedi, V. Pallottini, L. Canuti, M. De Canio, A. Urbani, V. Marzano, T. Cornetta, P. Stano, A. Giovanetti, L. Stella, A. Canini, G. Federici, and G. Ricci. Nuclear shield: a multi-enzyme task-force for nucleus protection. PLoS One, 5:e14125, 2010.

[32] L. Dock, M. Martinez, and B. Jernström. ncreased stability of (+/-)- 7 beta,8 alpha-dihydroxy-9 alpha,10 alpha-epoxy-7,8,9,10-tetrahydrobenzo [a]pyrene through interaction with subcellular fractions of rat liver. Chem.

Biol. Interact., 61:31–44, 1987.

[33] C.J. Roche, D. Zinger, N.E. Geacintov, and R.G. Harvey. Enhancement of stability of 7 beta,8 alpha-dihyroxy-9 alpha epoxybenzo(a)pyrene by complex formation with serum albumin. Cancer Biochem. Biophys., 8:35–40, 1985.

[34] S. Collins and M.A. Marletta. Purification of a benzo[a]pyrene binding protein by affinity chromatography and photoaffinity labeling. Biochemistry, 25:4322–4329, 1986.

[35] N. Merezhinskaya, G.A. Kuijpers, and Y. Raviv. Reversible penetration of alpha-glutathione S-transferase into biological membranes revealed by pho- tosensitized labelling in situ. Biochem. J., 335 (Pt 3):597–604, 1998.

[36] D. Mackay and S. Paterson. Evaluating the multimedia fate of organic chemicals: a level III fugacity model. Environ. Sci. Technol., 25:427–436, 1991.

[37] EPA. Soil screening Guidance: Technical background documents. Technical Report EPA/540/R95/128, US EPA, Washington, DC, 1996.

[38] J.M. Vanderkooi and J.B. Callis. Pyrene. A probe of lateral diffusion in the hydrophobic region of membranes. Biochem., 13(19):4000–4006, 1974.

[39] A.S. Verkman. Solute and macromolecule diffusion in cellular aqueous compartments. trends Biochem. Sci., 27(1):27–33, 2002.

[40] R.J. Ellis. Macromolecular crowding: An important but neglected aspect of the intracellular environment. Curr. Opin. Struct. Biol., 11(1):114–119, 2001.

(28)

[41] I.V. Tetko, J. Gasteiger, R. Todeschini, A. Mauri, D. Livingstone, P. Ertl, V.A. Palyulin, E.V. Radchenko, N.S. Zefirov, A.S. Makarenko, V.Y.

Tanchuk, and V.V. Prokopenko. Virtual computational chemistry laboratory - design and description. J. Comput. Aid. Mol. Des., 19:453–463, 2005.