Petri net–based model of the human DNA base excision repair pathway

Summary

Cellular DNA is daily exposed to several damaging agents causing a plethora of DNA lesions. As a first aid to restore DNA integrity, several enzymes got specialized in damage recognition and lesion removal during the process called base excision repair (BER). A large number of DNA damage types and several different readers of nucleic acids lesions during BER pathway as well as two sub-pathways were considered in the definition of a model using the Petri net framework. The intuitive graphical representation in combination with precise mathematical analysis methods are the strong advantages of the Petri net-based representation of biological processes and make Petri nets a promising approach for modeling and analysis of human BER. The reported results provide new information that will aid efforts to characterize in silico knockouts as well as help to predict the sensitivity of the cell with inactivated repair proteins to different types of DNA damage. The results can also help in identifying the by-passing pathways that may lead to lack of pronounced phenotypes associated with mutations in some of the proteins. This knowledge is very useful when DNA damage-inducing drugs are introduced for cancer therapy, and lack of DNA repair is desirable for tumor cell death.

Keywords:

DNA – Biology and life sciences – Cell biology – Chromosome biology – Chromatin – Chromatin modification – DNA methylation – Genetics – Epigenetics – DNA modification – Gene expression – Biochemistry – Nucleic acids – Computational biology – Research and analysis methods – Molecular biology – Molecular biology techniques – Molecular biology assays and analysis techniques – Biochemical simulations – DNA damage – DNA cleavage – DNA repair – Base excision repair – DNA synthesis – Ligation assay – Chemical synthesis – Biosynthetic techniques – Nucleic acid synthesis

Introduction

Living organisms as well as their functional parts such as organs, tissues, cells, etc. are very complex systems. These systems are composed of numerous basic building blocks connected by a dense network of interactions and mutual dependencies. A sum of such networks determines the structure and the functionality of the system. Therefore, in order to fully understand the nature of living organisms, it does not suffice to analyze properties of the building blocks, but it is also necessary to analyze the networks of interactions among them. Thus far, various computational methods for analyzing complex systems have been developed, and some of them can be applied in the analysis of biological systems [1, 2]. However, these systems have their own specificity, from which follows that in some cases it can be difficult to use methods developed earlier for other systems.

The first and necessary step in the analysis of a complex system is to build its formal model. This model can be expressed in a language of some branch of mathematics, and for instance, differential equations are often used for this purpose. Recently, Petri net theory has emerged as a promising approach for modeling and analysis of biological systems [1]. These nets have been described for the first time by C. A. Petri in the early 1960s in the context of theoretical computer science and for many years this was the main area of their applications [3, 4]. Since this time the theory of Petri nets has been intensively studied and many methods of the analysis of their properties have been developed. In the mid 1990s it was realized that Petri nets could also be used for modeling and analysis of biological systems [5, 6] and during the last two decades the possibilities and limitations of Petri net applications in this area were studied [1, 7, 8].

On the one hand, a strong advantage of Petri net-based models is their intuitive graphical representation. On the other hand, the nets can be analyzed using precise mathematical methods. Moreover, contrary to differential equations, they do not require the knowledge about precise values of some parameters of the modeled system, what can be a very important property in many cases since these values are often difficult to determine for biological phenomena.

The integrity of genetic information depends on the interplay between the DNA base modification (caused by enzymatic reactions), DNA damage (caused by cellular and environmental chemical reactions), and all DNA-binding proteins and readers of information that work together to make the whole system functional. The DNA damage repair is an exemplary complex biological system, in which a large number of elements (mostly enzymes that perform chemical reactions) cooperate with each other to protect the cell from consequences of DNA damage events. Base damage can be mutagenic and/or cytotoxic and can lead to mutagenesis or cell death, respectively, when unrepaired DNA lesions undergo replication or transcription. To avoid such consequences, correct repair of different types of DNA damage is assured by the existence of several DNA repair subsystems maintained by the cell in order to survive the DNA damage. In humans these subsystems, called pathways, include DNA damage signaling (DDS), direct reversal repair (DRR), base excision repair (BER), nucleotide excision repair (NER), mismatch repair (MMR), homologous recombination repair (HRR), nonhomologous end-joining (NHEJ) and translesion synthesis (TLS). Details concerning DNA damage and repair in several living organisms were collected and included in the public databases developed by our groups (http://repairtoire.genesilico.pl/pathways/ [9], https://www.ibb.waw.pl/pl/Sekretariat%20Naukowy/Projekty%20strukturalne/dnatraffic [10]).

The DNA repair processes have already been a subject of computational modeling. A broad range of modeling approaches gave opportunities to gain a variety of insights into DNA damage and repair. Some of these models were built based on detailed kinetic data available from experimental studies. Examples include the stochastic model of NHEJ, which uses molecule numbers, reaction rate constants and kinetic rate constants to explain the dynamic response to damage induced by different levels of gamma irradiation in human fibroblasts [11]. Computational modeling was also conducted for some sub-pathways of BER, where e.g., differential equations were used to model the removal of 8-oxoguanine (8-oxoG) via base excision repair initiated by human OGG1 glycosylase [12]. In this case the availability of kinetic data for both wild type enzymes and their mutant variants allowed for the prediction of the kinetics and capacity of the whole sub-pathway in the presence of some single-nucleotide polymorphisms which have been associated with cancer. Other examples of kinetic models for BER include the works of [13] and [14]. Valuable insights into DNA repair processes were also delivered by models which did not require incorporation of very detailed data as in the cases discussed above. For example, Monte Carlo simulations of radiation-induced damage and its repair were used to investigate the repair of clustered damage via BER and NER pathways [15, 16]. These approaches allowed to estimate the probabilities of correct repair, fractions of damages converted into double-strand breaks (DSBs) and fractions of repairs ending with mutations.

In this work, we present an extensive study of a Petri net model of the BER repair pathways, which we have built based on available literature and from biological pathway databases. In contrast to models built by others, our approach involves modeling the repair of many different types of DNA damage through the activity of a collection of DNA glycosylases initiating BER and other enzymes performing further steps of the repair. With such a comprehensive model formal analysis becomes possible using methods provided by the Petri nets theory. Here we have performed a t-invariant based analysis, involving MCT sets and t-clusters, decomposing the model into precisely identified and named functional modules. Understanding their interaction provides valuable insights into the complex DNA repair process. We have also performed an in silico knockout analysis as a prediction of system behavior under a collection of disturbances.

The organization of the paper is as follows. Section 2 provides an overview of Petri nets, including the analytical analysis of t-invariants. In Section 3 the results of the net analysis are presented. The paper ends with discussion and conclusions given in Section 4.

Methods

Petri nets

A Petri net is a mathematical object based on a weighted directed bipartite graph. There are two disjoint sets of vertices called places and transitions. Places correspond to passive components of a modeled system (e.g., chemical compounds) while transitions model active components (e.g., chemical reactions or other elementary subprocesses). Places can be connected with transitions and transitions with places by arcs, which describe causal relations between elementary components of the system. The arcs are labeled by weights being positive integer numbers. Places can hold objects called tokens, which represent discrete quantities of the passive components. In a graphical representation of Petri nets places are denoted as circles, transitions as rectangles or bars, arcs as arrows and tokens as dots or positive integer numbers within places. If a weight of an arc is equal to one, usually it is not shown in the graphical representation; otherwise it is represented as a number labeling the arc.

The bipartite graph underlying a Petri net determines its structure (which should correspond to a structure of the modeled system), but most of the fundamental properties of nets of this type follow from their dynamics which is related to tokens. The flow of tokens through the net from one place to another via transitions corresponds to the flow of information, substances etc. through the modeled system. The distribution of tokens over all places, called marking, represents a state of the system.

The flow of tokens is governed by a transition firing rule. According to this rule, a transition is enabled if the number of tokens in each of its pre-places, i.e., the ones which directly precede this transition, is equal to at least the weight of the arc connecting such a place with the transition. An enabled transition may fire, which means that tokens flow from its pre-places to its post-places, i.e., the ones which are its immediate successors. The number of flowing tokens is equal to the weight of the arc. There are two exceptions to this rule, i.e., transitions without pre-places are continuously enabled while transitions without post-places do not produce tokens. These types of transitions can be used to model interactions of the system with its environment. A read arc is a special kind of bidirectional arc, corresponding to a pairs of arcs. Tokens in a place connected by a read arc to a transition are not consumed when firing this transition, but they are required to enable it.

Part of the results in the following paper comes from the analysis of the net simulation. In such a simulation one must set specific rules governing the firing of the transitions. Our model is based on a classical Petri net, where in general only three scenarios for transition firing are available: in a single step of the simulation only one of the enabled transitions can fire, all enabled transitions can fire simultaneously or, in a third case, only some, randomly chosen enabled transitions will fire. The latter scenario has been used in our simulations. Specifically, each enabled transition has 50% chance of firing and the sequence of enabled transitions selected for firing in a single simulation step is (each time) chosen randomly.

More complex firing scenarios could be available if, e.g., stochastic Petri net has been used. It would however require assigning a so called firing rates to each transition, what in general would be a difficult task requiring knowledge about more or less specific probability of each reaction involved in the modeled system. The assumption of 50% firing chance is therefore a trade-off, allowing simple yet still valuable stochastic simulations of the model.

As a summary one can say that in a single step some random enabled transitions (with 50% chance for firing) are randomly assigned to the firing list. One must note that they will not necessarily all fire, because competing transitions on such a list can share the same pre-places and consume “activation” tokens from each other depending on the firing sequence (i.e., when a transition is on that list, but because of firing of the other ones such a transition loses the necessary minimum of tokens in its pre-places, it stops being active and cannot fire). However, such a firing sequence is random in each simulation step, so a sufficient stochastic scenario is provided to study the dynamic behavior of the model.

More advanced analysis methods for simulation of a net would be available if a stochastic Petri net (SPN) or continuous Petri net (CPN) have been used to model the DNA repair process. The latter would require the formulation of ODEs for the reaction flows. For the SPN model, the so-called firing rates need to be assigned for each transition—in the stochastic Petri net they influence the transition chances of firing in the simulation. Both approaches have advantages (e.g., a more complex and accurate simulation) and disadvantages. As for the latter, it is still difficult to find specific and accurate values which could describe the ODE in the CPN model or firing rates for all transitions for the SPN model. Further, in a continuous Petri net an invariant, MCT and cluster-based analysis explained in the next section is not available, while in this work such an analysis allowed us to divide the Petri net model into some functional subnets, characterized by the common functionality.

t-invariants analysis

Besides the graphical representation of Petri nets also a more formal one, called incidence matrix, is used. Entries of such matrix A = [a_ij]_{n × m} are integer numbers, and entry a_ij is equal to a difference between numbers of tokens in places p_i before and after firing transition t_j.

Among many properties of Petri nets, the ones related to t-invariants are especially important in the context of an analysis of biological systems models. A t-invariant is a vector x of integer numbers being a solution to the equation

t-invariants correspond to some subprocesses of the modeled biological system which do not change its state. Hence, analyzing relations among them may lead to discoveries of properties of the system [17]. The net is covered by t-invariants if every transition belongs to a support of at least one of them. In such a situation every transition contributes to at least one process of the modeled biological system [18–20]. Therefore, when an analysis is based on t-invariants the net should be covered by them. In order to find relations between t-invariants they are grouped into sets called t-clusters, using standard clustering algorithms. Also transitions can be grouped into sets called maximal common transition sets (MCT sets) what helps in the analysis of t-invariants [17, 21]. Both t-clusters and MCT sets correspond to some functional modules of the biological system. A more precise description of a Petri net elements, invariants and MCT sets is available as Appendix in S1 Appendix.

Results

The model of the human base excision repair pathway

The Petri net presented in this work is a model of the human base excision repair pathway. The model was built based on data from the publicly available databases: REPAIRtoire—a database of DNA repair pathways [9], DNAtraffic [10] and Reactome [22, 23]. Subsequently, the computational model was supplemented with several new aspects of the apurinic/apyrimidinic (AP) site processing which are not included in the above-mentioned databases but were reported in the literature:

• displacement of OGG1 glycosylase at the AP site by AP endonuclease [24–26] or NEIL1 glycosylase [27, 28],

• displacement of NTH1 glycosylase by AP endonuclease [29],

• Polβ-mediated LP-BER [30, 31],

• mechanism of the repair initiated by NEIL2 and NEIL3 [32, 33].

The model has been built and analyzed using Holmes software [34] and is available in SPPED (Snoopy), XML (SMBL), .project (Holmes) and PDF formats as S1, S2, S3 and S4 Files respectively.

The model can be subdivided into 4 subnets, described below:

Subnet 1: Introduction of DNA damage

This part of the model represents the creation of different types of DNA damage from undamaged DNA.

Subnet 2: AP site formation

The biggest part of the model represents the recognition of DNA damage by DNA glycosylases, followed by the removal of damaged bases by cleavage. This part of the model can be further subdivided into 11 modules, each of which represents the activity of one DNA glycosylase (Table 1).

DNA glycosylases and their substrates are represented by places. The substrates were divided into groups, based on a subset of DNA glycosylases that recognize them (Table 2).

Each group of DNA damage is represented by a separate place. For each glycosylase-substrate pair this part of the model contains a transition representing the damage recognition process, which leads to the glycosylase-damaged DNA complex formation and cleavage of the damaged base, which in turn results in the AP site formation. Fig 1 presents the submodule for the OGG1 DNA glycosylase (full names of places and transitions of this submodule are given in Table 3).

Subnet 3: AP site processing

AP site (generated by damaged base removal) is incised either by AP-lyase activity of a bifunctional DNA glycosylase or by AP endonuclease (APE). AP-lyase activities of OGG1, NTH1, NEIL1, NEIL2 and NEIL3 are represented in the model by appropriate lyase activity transitions. In the case of OGG1, NTH1 and NEIL3, which perform β-elimination, the lyase activity results in the formation of DNA incision with 3ʹdRP (deoxyribose phosphate) end (Fig 2A). Full names of places and transitions for all three Fig 2 subnets are given in Table 4.

Only NEIL1 and NEIL2 generate DNA incision with 3ʹP end through β,δ-elimination (Fig 2B). In contrast, AP sites generated by monofunctional glycosylases (SMUG, MBD4, MPG, MYH, TDG, UDG) need additional processing and involvement of APE. This reaction results in DNA incision with 5ʹdRP end (Fig 2C).

The AP-lyase activity of OGG1 is weak ‒ after base cleavage, OGG1 can be substituted by APE, which catalyzes strand incision leading to 5ʹdRP end (Fig 2A). Similarly, APE can process AP sites generated by NTH1, and NEIL1 can catalyze β,δ-elimination at the AP sites generated by OGG1 (Fig 2A and 2C).

The 3ʹdRP and 3ʹP DNA termini (dirty ends) are refractory to DNA polymerase repair synthesis and need to be “cleaned” to 3ʹOH for BER to proceed. In human cells, APE and PNKP are responsible for the removal of 3ʹdRP and 3ʹP, respectively. 5ʹdRP end prevents the ligation process and its removal is primarily carried by DNA Polβ or FEN1.

Subnet 4: DNA repair synthesis and ligation

The last step of BER is the replacement of the excised nucleoside by repair synthesis catalyzed by DNA polymerases (presented in Fig 3, full names of places and transitions are given in Table 5) and nick sealing by a DNA ligase. It can proceed either as SP- or LP-BER. Typically, SP-BER is performed by Polβ and LIG3-XRCC1 complex on the products of bifunctional glycosylases activity. Repair initiated by monofunctional glycosylases proceeds via LP-BER which requires activity and presence of Polβ, Polδ/Polε, PCNA, FEN1 and LIG1. However, the monofunctional glycosylase pathway can also proceed via SP-BER, after removal of 5ʹdRP by Polβ or FEN1. If Polβ removes 5ʹdRP than the final ligation step is performed by LIG3-XRCC1 complex, otherwise it is done by LIG1. DNA synthesis in LP-BER is usually performed either by Polδ or Polε, with Polδ being preferred in high concentrations of PCNA [66]. However, it was also shown that Polβ can conduct strand synthesis in LP-BER.

Mathematical analysis of the model

The presented Petri net model consists of 259 transitions and 179 places. Full names of places and transitions are given in Tables A and B in S1 Tables, respectively.

In the first step of the analysis of the model, some of its structural properties have been identified. The net is pure and ordinary, i.e., there are no read arcs, and all the arcs have a weight equal to 1. Therefore, the net is homogenous–all outgoing arcs of a given place have the same weight. The network is connected, but not strongly connected, meaning that all its nodes have non-directed connection(s) with the others, but there are pairs of nodes without a directed path between them. The net is not conservative, i.e., many transitions consume a different number of tokens from their pre-places than they produce in the post-places. There are also static conflicts–some transitions share the same pre-places. Finally, there are input and output transitions in the model, which correspond to the connections of the analyzed system with other processes.

An important part of the analysis of the net is based on the t-invariants. There are 245 t-invariants and the net is covered by them. In the presented model of DNA BER-repair process, there are 94 non-trivial MCT sets, i.e., containing more than one transition. Many of the MCT sets are responsible for almost the same functionality, but operating on different substrates. Description of the most important MCT sets is presented in Table 6.

The most important basic functions of the modeled system correspond to eight MCT sets: m₁ (short patch repair with Polβ and PNKP), m₁₄ followed by m₁₅ (short patch repair for DNA-APE complex), m₁₆, m₁₇, m₄₀ (short patch (SP) and long patch (LP) repair in which Polβ and FEN1 are involved), finally m₂ and m₃ responsible for long patch repair using Polδ and Polε, respectively. Two very important sets are m₁₃ and m₁₄. The latter is crucial for starting one of the short patch repairs with Polβ. Set m₁₃ is responsible for the preparation of the DNA complex for four out of six repair paths, including three long patch repairs (with Polβ, Polδ, and Polε). Overall, the results of analysis of the MCT sets suggest that in most cases the chosen repair sub-pathway is determined by the damage type at the very beginning of the process. However, in some cases alternative pathways are available. For example repair processes in which generation of 5′dRP by APE activity takes place, the repair can either proceed via SP-BER (m15) or via LP-BER, if PCNA level is high (m17, m40).

The next step of the model analysis concerns t-clusters computation. Such successfully obtained clusters can divide the Petri net into different subnets responsible for the similar or common functions in the whole DNA repair process. Such a division is based on the t-invariants set. The description of such subnets/clusters can use the previously computed MCT sets for simplification (by replacing some groups of transitions by such sets), yet one should have in mind that such a t-cluster analysis is based solely on the computed t-invariants sets.

In order to find the optimal number of clusters, many similarity measures and grouping algorithms had to be used and their results compared. The theory is that for each cluster (group) its t-invariants must be very similar to each other according to the chosen similarity measure, and as much as possible dissimilar to t-invariants belonging to different groups [67]. Three important choices needed to be made: choosing the proper similarity measure, the clustering algorithms and the number of clusters. Thorough comparative tests of different choices in these three variables have been performed. In the tests eight similarity measures have been used: Binary, Canberra, Euclidean, Manhattan, Maximum, Minkowski, uncentered and centered Person (the latter often called Correlated Pearson metric). Further, seven clustering algorithms have been tested: Centroid, Complete, McQuitty, Median, Single, UPGMA, and Ward. The evaluation concerned the number of clusters ranging from 2 to 30. In order to evaluate the results two indices have been used, i.e., Mean Split Silhouette (MSS) [68] and Caliński-Harabasz index [69]. For the presented Petri net the best clustering has been obtained using the UPGMA algorithm and correlated Pearson similarity measure. The best number of clusters according to the used evaluation measures is 27. However, 20 of them are trivial single t-invariants responsible for the fluctuations of the system components, like APE, Polβ, Polδ, Polε and the glycosylases responsible for damage detection. The obtained clustering (i.e., set of clusters) is presented in Table 7.

Clusters 21–27 from Table 7 are all responsible for different repair paths within the presented net. Table 8 presents a more detailed description of each path, providing information about e.g., glycosylases involved in damage detection or main polymerases involved in the repair process.

Table 9 contains the names of the most important transitions corresponding to the elementary processes within short and long patch repair pathways.

Analysis of the robustness of the BER pathway via systematic in silico knockout experiment

The impact of each modeled repair synthesis sub-paths: three short-patch (SP) ones (where DNA Polβ is always present) and three long-patch (LP) repairs (with polymerases δ, ε or β) were analyzed via systematic in silico knockout experiment. The six repair synthesis sub-paths defined for this analysis are as follows:

• Path 1: SP₍₁₎, refers to t-invariants cluster 21, DNA synthesis performed by Polβ, includes PNKP,

• Path 2: SP₍₂₎, refers to t-invariants cluster 22, DNA synthesis performed by Polβ,

• Path 3: LP₍₁₎, refers to t-invariants cluster 23, DNA synthesis performed by Polβ,

• Path 4: SP₍₄₎, refers to t-invariants cluster 27, DNA synthesis performed by Polβ, includes FEN1,

• Path 5: LP₍₂₎, refers to t-invariants cluster 25, DNA synthesis performed by Polδ,

• Path 6: LP₍₃₎, refers to t-invariants cluster 26, DNA synthesis performed by Polε.

Net simulations were used to perform a series of knockout experiments. Each simulation contained 20,000 steps and was repeated 20 times. In a single step, any enabled transition had a 50% chance for firing. Transition t₁₆₉ (DNA_back_to_pool) was treated as a marker of successful DNA repair. Its chance of firing represents the chance of successful completion of the repair process. The results of the simulations for the non-disturbed model are summarized in Table 10. The chance of successful repair, represented by the firing of transition t₁₆₉, equals on average 23.3%. The biggest proportion of successful repairs is processed via the SP₍₂₎ pathway.

The results of the knockout experiments are presented in Table 11. Transitions critical for each repair synthesis sub-paths have been identified, and their impact on the overall chance of successful repair was checked. The most prominent change in the repair success rate was observed when Path 2 was disabled. In this setup Paths 1 and 3–6 exhibit a slightly higher activity but they cannot fully recapitulate the repair. At the same time the lack of activity in Paths 3–6 did not result in an observable change in DNA repair efficiency. The high importance of Polβ in DNA repair is a well-known phenomenon—the lower repair rate is thus expected (see Discussion for more details). To the best of our knowledge, Polδ and Polε knockouts influence on BER efficiency and higher activity of other BER sub-pathways after Polβ knockout have not been reported yet.

We have also investigated influence of a glycosylase knockout on the modeled behavior of BER, taking the TDG glycosylase as an example. We found that TDG knockout, simulated by inactivation of transition t_206, which introduces TDG into the system, leads to: i) higher activity of glycosylases involved in processing of the same DNA damage substrates as TDG (namely MBD4, OGG1, NEIL1, NEIL2 and SMUG), ii) slight decrease in the overall DNA repair efficiency, iii) accumulation of DNA damage from group 11, which are recognized only by TDG and iv) decreased activity of the SP-BER pathway.

Discussion

The BER process can be described as a unique DNA repair pathway, where very precise enzymes are used to recognize and remove specific DNA damage. Thus far eleven DNA glycosylases have been found in human cells (Table 1). All DNA glycosylases are divided into four structurally different superfamilies and share a similar ability to recognize the DNA damage. However, there are also several mismatches (e.g., T-C) or DNA damage types (e.g., 3-mA, 7-mG) which are removed only by one human enzyme (Table 2) [9, 10, 36, 65, 70].

Although BER has already been a subject of several computational modeling attempts, we provide the most comprehensive model, which includes the widest range of protein factors involved in the DNA repair process and the biggest collection of DNA damage types (Table 1). The models available to date concentrate mainly on the repair initiated by the OGG1 glycosylase [12–14], while our computational model of human BER includes 68 different DNA damage types found in living cells and identified and repaired by 11 human DNA glycosylases. We have also included into our Petri net model several mechanisms that have not been considered in the available models, e.g., displacement of OGG1 at the AP site by AP endonuclease or NEIL1, displacement of NTH1 by AP endonuclease, Polβ-mediated LP-BER and the mechanism of the repair initiated by NEIL2 and NEIL3. However, we have included in our model mainly the major DNA damage types. As a consequence, some of the products of minor G base oxidative damage (e.g. 5-carboxamido-5-formamido-2-iminohydantoin, 2Ih) which can be removed from DNA by hNEIL1-3 glycosylases [71] are not present in the model. We have also omitted several bulky DNA damage products (from the minor fraction) which are used mostly for in vitro studies of DNA glycosylases activities.

BER can proceed via different sub-pathways, depending on the damage and the type of DNA glycosylase that detected the problem (Table 2). In general, the repair DNA synthesis can be performed by DNA Polβ (Polβ, used for both the so-called short-patch (SP) and long-patch (LP) base excision repairs). Two other major DNA polymerases δ and ε (Polδ and Polε) can repair DNA instead of Polβ, and they are always involved in the LP-BER pathway with PCNA involvement. Other possibilities for repairing DNA by the BER mechanism have also been included in our model. The analysis performed on the base of net t-invariants allowed to divide a structure of the net into meaningful biological units (MCT sets) and more importantly to group t-invariants which represent basic subprocesses of the model into clusters (Table 6). The clusters allow for separating different repairing processes, showing the glycosylases, polymerases and other compounds involved in the repair processes, as presented in Tables 7 and 8 and Figs 1–3.

Using in silico simulations, we showed that most efficient and successful DNA repair of the studied DNA damage was reached via SP-BER (Table 10). This result was consistent with biological data showing that the patch size in the nucleosomal BER was rather shorter (1 nt) than longer (2–12 nt) [72]. Our in silico knockout experiments showed that Polβ, the so-called repair polymerase (Table 11, Path 2) is necessary for BER. Lack of DNA Polβ in our model resulted in stopping the cleaning of 5’dRP moiety and gap filling (DNA synthesis) which was simultaneously manifested as a block of whole SP-BER. In living cells, the perturbation in BER can lead to genomic instability. Our results could be supported by the observation that somatic mutations in Polβ gene have been detected in various types of cancers (NCI Genomic Data Commons (GDC)). However, it is still not clear whether and how Polβ mutations and overexpression can be linked to cancer onset and its progression. Recently, it has been shown that human R152C Polβ mutant was impaired in BER activity and efficiency. Moreover, the mutant cells have displayed a high frequency of chromatid breakages [73]. Mice carrying a targeted disruption of the Polβ gene has shown growth retardation and died of a respiratory failure immediately after the birth [74]. Also, homozygous Polβ R137Q knock-in mice embryos were typically small in size and had a high mortality rate (21%). In this mutant the BER efficiency was impaired, which subsequently ended in double-strand breaks (DSBs) and chromosomal aberrations [75]. Recently it has been shown that a mouse model with decreased expression of Polβ (with Y265C mutation) developed systemic lupus erythematosus (SLE) [76].

It is not clear which of the two activities of Polβ is the most important in SP-BER. Both DNA synthesis and 5′-dRP-lyase activity can be replaced by other DNA polymerases or FEN1, respectively. If the 5′-dRP moiety is reduced or oxidized, the 5′-dRP-lyase of Polβ cannot remove the modified sugar residue, and LP-BER is initiated. Next, a flap (2–10 nucleotide long) is subsequently removed by FEN1. FEN1 mutations are very rare, suggesting that FEN1 is important for normal DNA metabolism [77, 78]. Targeted deletion of the Fen1 gene in mice causes early embryonic lethality [79]. However, several somatic FEN1 mutations have been detected in human cancer (GDC) but the relationship between FEN1 deficiency and cancer susceptibility remains unclear. Recently, it was reported that L209P FEN1 mutation is associated with colorectal cancer. Human L209P FEN1 mutant was lacking the exo- and endonuclease activities but retaining DNA-binding affinity. This was a dominant-negative mutation and mutated protein impaired LP-BER in vitro and in vivo [80]. In contrast, our computational Petri net BER-knockout model with eliminated FEN1 activity did not show any significant change in BER efficiency (Table 11, Path 4). This inconsistency can be explained by the fact that FEN1 functions not only in LP-BER and it is mostly recognized as a central component of cellular DNA metabolism (e.g. processing of Okazaki fragment maturation intermediates, telomere maintenance and rescue of stalled replication fork) [77].

It is estimated that each day as many as 10,000 abasic sites are formed in one human cell [81]. Removing AP-sites from DNA is a daily task for the DNA repair/tolerance system. Apurinic/apyrimidinic endonuclease 1 is responsible for the AP sites processing which is necessary for further steps of DNA repair pathways [82, 83]. In our computational BER model, the elimination of AP-sites endonuclease 1 (APE1) resulted in downstream blocking of DNA repair, both, SP-BER and LP-BER. Lack of AP-sites repair can be an explanation of the embryonic lethality of Ape1 null mice [84, 85] and increasing tumor susceptibility in heterozygous mice [86]. APE1 mutations in humans (APE1 variants: L104R, E126D, and R237A, exhibiting approximately 40–60% reduction in specific incision activity) have been associated with amyotrophic lateral sclerosis (ALS) [87, 88]. Also, somatic APE1 mutations (APEX P112L, W188X, and R237C) were found in endometrial cancers [89].

Preparation of the 3′OH ends is a key step for repair DNA synthesis. In our model, we observed a reduction of a chance for successful repair when polynucleotide kinase/phosphatase (PNKP) was removed from the pathway (Table 11, Path 1). PNKP involvement in BER is mostly due to APE-independent base excision repair pathway in human cells after NEIL1 and NEIL2 action on oxidized base lesions [27, 28]. However, PNKP (two activities: DNA 5'-kinase and DNA 3′-phosphatase) is mostly recognized as an enzyme which generates 5′-phosphate/3′-hydroxyl DNA termini that are critical for ligation by the non-homologous end joining (NHEJ) DNA ligase LigIV during double-strand break repair (DSBR). Microcephaly with early-onset, intractable seizures and developmental delay (MCSZ) is a hereditary disease caused by mutations in PNKP [90, 91].

Eleven human DNA glycosylases were in silico deleted one-by-one in our computational BER model. None of them was critical for BER function because almost each DNA damage can be removed by at least two distinct glycosylases (Table 2). Apart from DNA damage removal, various BER glycosylases (SMUG, MBD4, TDG, UNG2, NEIL1-3) are involved in nucleotide replacement during the active DNA demethylation process [57]. However, unlike other DNA glycosylases, TDG is essential for embryonic development; mice die at day 11.5 [92]. The lethal phenotype is associated with epigenetic aberrations affecting the expression of developmental genes. Mouse embryonic fibroblasts, (MEFs) derived from Tdg null embryos showed impaired gene regulation due to imbalance histone modification and CpG methylation [93]. Other DNA glycosylases are not essential, mouse knockouts are viable, however, showing increase in mutation frequency or some immune dysfunction [36]. Human mutant DNA glycosylases have also been widely described and correlated with a predisposition to various diseases [70]. Also, an interplay between DNA repair of the oxidatively damaged base, 8-oxoG (8-oxo-7,8-dihydroguanine) and transcriptional activation has been documented for mammalian genes: removal of 8-oxoG from the coding strand by OGG1-mediated BER resulted in upregulated transcription [94]. Since oxidation is ongoing and transition of C to U occurs spontaneously or at specific times during differentiation and development, there is a strong suggestion that BER substrates might be epigenetic and modulate transcription factor binding [71, 95].

We believe that our model can be a valuable tool that can help to predict the influence of changes in protein activities or levels on the repair process as well as some regulatory functions. It can be used to predict the sensitivity of the cell with inactivated repair proteins to different types of DNA damage and it can help to identify the by-passing pathways that may lead to lack of pronounced phenotypes associated with mutations in some of the proteins.

Supporting information

S1 Tables [pdf]

S1 Appendix [pdf]
Appendix with more detailed description of Petri net elements and theory.

S1 File [spped]
The model in SPPED format.

S2 File [xml]
The model in SMBL format.

S3 File []
The model as Holmes project.

S4 File [eps]
The model as a figure in PDF format.