DNA molecules packaged in our chromosomes carry our genetic blueprint, and their preservation is essential for the coordination of cellular function and organization of life (Branzei and Foiani, 2008). As DNA is the repository of our genetic information, we would expect its structure to be highly stable. This is not the case for the DNA (Reed and Waters, 2005). Damage to DNA is a constant threat (Lindahl, 1993; Alberts, 2015). The DNA molecule is intrinsically reactive, as it is very susceptible to chemical and physical factors, which can lead to DNA lesions, such as base loss, base modification, and double-strand DNA breaks (Hoeijmakers, 2009; Çağlayan and Wilson, 2015; Ross and Truant, 2016; Yadav et al., 2020). Physiological conditions such as oxygen-rich, aqueous, or pH 7.4 (Lindahl, 1993) as well as chemical events such as hydrolysis and exposure to reactive oxygen species (ROS) or other reactive metabolites can damage DNA. Exogenous chemicals or endogenous metabolic processes trigger chemical reactions. Although exogenous stressors can be extremely powerful, endogenous threats are constant and unabating. It is estimated that a single cell experiences up to 10⁵ spontaneous or induced DNA lesions per day (Lindahl, 1993; Bont, 2004; Kovalchuk, 2016; Chatterjee and Walker, 2017).

DNA damage has far-reaching consequences, such as preventing RNA polymerase from transcribing the correct messenger RNA sequence to produce the correct protein. In the longer term, cellular malfunctions such as cancer initiation, inborn defects, and ageing that result after damaged DNA replicates are examples of unpredictable long-term consequences of DNA damage; as base misincorporation causes mutations which alter the genetic code (Reed and Waters, 2005). Therefore, a coordinated response to DNA damage is necessary in order to ensure cellular viability and prevent diseases. Cells, fortunately, possess a robust system of mechanisms that function together to reduce the adverse consequences of DNA damage and ensure that their genetic information is faithfully replicated, thus maintaining the integrity of their genome (Ganai and Johansson, 2016). This coordinated effort, known as DNA damage response (DDR) operates by sensing and signalling the genotoxic events, and the damage is then resolved either by DNA repair machineries, or cell death if DNA cannot be repaired. DNA repair functions as part of the DNA damage response (DDR) (Liu, 2001; Chatterjee and Walker, 2017; Reed and Waters, 2005; Hoeijmakers, 2001, 2009).

DNA repair has so far been shown to exist in both prokaryotic and eukaryotic organisms, with over 150 proteins directly involved in safeguarding the genome (Sancar et al., 2004; Friedberg et al., 2005; Wood et al., 2005; Yousefzadeh et al., 2021). DNA repair processes restore DNA back to its normal sequence and structure after damage (Friedberg et al., 2006), and are characterised traditionally by the type of damage they repair. There are five major DNA repair pathways available to cells to deal with DNA damage burdens. Each of these processes recognises a particular type of DNA lesions, and together work in preventing mutagenesis. They include 1) direct reversal repairs that repairs lesion induced mainly by alkylating agents, 2) Base excision repair (BER), for small base modifications like single-strand breaks (SSBs) and non-bulky damaged DNA bases, 3) Nucleotide excision repair (NER), that corrects bulky, helix-distorting DNA lesions, 4) mismatch repair (MMR), that repairs base-base mismatch and insertion or deletion loops (IDLs), 5) Recombinational repair, which is divided into non-homologous end joining (NHEJ) and homologous recombination repair (HR), both of which repairs DNA double-strand breaks (DSBs). Other types of DSB repairs include alternative non-homologous end-joining (alt-NHEJ, MMEJ) and translesion synthesis (TLS), which operates as a tolerance mechanism for DNA damage (Jackson and Bartek, 2009; Hosoya and Miyagawa, 2014; Li et al., 2021).

For frequently occurring DNA damage, direct reversal of DNA damage by specialised proteins is the most efficient and most straightforward method of DNA repair. However, this approach is only used by a small proportion of DNA repair types. Most damage to DNA is repaired by the removal of damaged bases and is followed by resynthesis of the removed/excised region (replacement) (Cooper, 2000). The pathways involved in the removal of base damage are base excision repair (BER), nucleotide excision repair (NER) and mismatch repair (MMR) (Cooper, 2000). The rest of the pathways repair damage to DNA structure/backbone. DNA damage can cause breaks in the DNA backbone, single-strand breaks in one strand, or double-strand breaks on both strands. Single-strand breaks are repaired by mechanisms sharing common steps in the BER pathway; however, DSBs are especially harmful as, by definition, no unbroken complementary strand exit which can serve as a template for repair when both strands break (Bennett et al., 1993; Friedberg et al., 2006). For cells with DNA already replicated prior to cell division, the duplicate copy can easily supply the missing information. So in these cells, DSBs can be repaired by HR, involving the exchange of DNA strands.

Even so, very efficient repair mechanisms can sometimes fail to provide a clean template for DNA synthesis. Replication errors can make it past these mechanisms, as DNA repair can also undergo mutations and become dysregulated. DNA repair gene mutations have been known to cause a variety of rare inherited human syndromes. Some of which include premature ageing phenotypes, increased sensitivity to ionising radiation exposure, and increased cancer risk (Friedberg et al., 2006; Lok and Powell, 2012; Carusillo and Mussolino, 2020). It has also been found that inherited defects in each of the DNA repair pathways are associated with distinct genome instability syndromes (Yousefzadeh et al., 2021), syndromes characterised by developmental defects (Bouwman and Jonkers, 2012; Ghosal and Chen, 2013; Wolters and Schumacher, 2013; Wood, 2018).

The dysregulation of DNA repair gene networks underlies many human genetic diseases that affect a wide range of body systems but all share a common trait, predisposition to cancer (Chatterjee and Walker, 2017). Almost all human cancers are spontaneous, not inherited, and are caused by environmental or genetic factors. It is of great public health interest to determine which genetic variations increase cancer risk in normal populations, and DNA repair genes are likely contenders. Therefore, elucidating the molecular mechanism behind DNA repair defects may provide a framework for understanding the complex pattern of genetic variations that contribute to spontaneous human cancers.

In this study, we demonstrate how to translate a network (mathematically definable as a hyper-graph) into a set of first-order differential equations of the mass action law type. Once the model has been established, we present its numerical solution and carry out a sensitivity analysis of its kinetic rates, whose numerical values are mostly unknown. The results of the analysis of network dynamics complement those produced by the calculation of centrality measures, and together they produce a set of genes of similar biological interest, and in perspective also of medical interest, due to their characteristics of topological centrality and vulnerability to stimuli. The paper is organized as follows: in Section 2 we describe the mechanisms of double-strand break repair pathway homologous recombination repair necessary to understand and interpret the results of the computational analysis, in Section 3 we describe the rules on which the automatic translation of the network into a rate equation system is based and the methods of sensitivity analysis of the model. In Section 4 we present the results of the analysis, and, finally, in Section 5, we draw some conclusions.

DSBs are the most serious DNA damage, as both DNA strands are impaired simultaneously. Therefore, due to the magnitude of differing factors leading to DSBs, the effectiveness of their repair is crucial for cell survival and the functioning and prevention of DNA fragmentation, chromosomal translocation and deletion. DSBs can be repaired in mammalian cells by NHEJ, HR, and single-strand annealing (SSA). Unrepaired SSBs result in much more cytotoxic DSBs formation during the S-phase progression of the cell cycle (Kennedy and D’Andrea, 2006). Homologous recombination is a process by which DSBs are repaired through the alignment of homologous sequences of DNA (Dietlein and Reinhardt, 2014) and occurs primarily during the late S to G2 phase of the cell cycle (Cerbinskaite et al., 2012; Chatterjee and Walker, 2017).

Homologous recombination is the second major DSB repair pathway and requires a second, homologous DNA sequence to function as donor template. There are two phases to this process, the first phase triggered by sensor proteins that belong to the MRN complex, and the second phase by the stimulation of resection steps, initiated in the first phase and subsequently extended. HR generally involves the following stages:

1. DSBs are recognised and sensed by the MRN complex (Kim et al., 1994), which activates ATM kinase, initiating the DSB end resection steps, where CtBP-interacting protein (CtIP) and the MRN complex work together to generate single-strand DNA (ssDNA) at the DSB ends ((Zhao et al., 2020).

Following the recognition and sensing of DSBs, a process known as DNA end resection is activated, a critical function in HR (Liu and Huang, 2016; Zhao et al., 2020). DNA end resection catalyses the nucleolytic degradation of the broken ends of DSBs (by the CtIPMRN complex) in the 5′ to 3′ direction generating 3′ single-stranded DNA (ssDNA). The 3′ ssDNA then provides a platform for the recruitment of HR repair-related proteins (Huertas, 2010; Liu and Huang, 2016; Zhao et al., 2020). Following the generation of ssDNA, downstream nucleases and helicases, such as exonuclease 1 (EXO1) or DNA replication ATP-dependent helicase/nuclease DNA replication helicasenuclease 2 (DNA2) and Bloom syndrome protein (BLM), are conscripted to extend the 3’ ssDNA for HR repair (Huertas and Jackson, 2009; Yun and Hiom, 2009; Zhao et al., 2020). The identities of these DNA helicases and nucleases are yet to be clearly defined in humans (as in yeast), partly because there are many candidate proteins. Although five RecQ helicase homologs have so far been identified in yeast (Bloom helicase [BLM], Werner helicase/nuclease [WRN], RECQ1, RECQ4, and RECQ5) (Chu and Hickson, 2009; Lu and Davis, 2021), convincing evidence point up BLM in resection (Gravel et al., 2008; Nimonkar et al., 2008). Following resection, the exposed single-strand DNA (ssDNA) is recognised and bound by RPA complex for protection.

RPA plays a significant role in coordinating DNA resection processes and simultaneously preserving the integrity of the resultant ssDNA (Sun et al., 2019). RPA is a heterotrimeric ssDNA binding protein essential to nearly all DNA processing events and associates with ssDNA with very high affinity (Kd sim 109–1010 M) (Maréchal and Zou, 2014). It is comprised of three protein subunits, RPA70, RPA32 and RPA14 and contains multiple oligonucleotideoligosaccharide (OB)-folds that interact with both ssDNA and proteins (Kim et al., 1994; Fanning, 2006; Feldkamp et al., 2014; Maréchal and Zou, 2014). RPA is flexible (Brosey et al., 2013). Its versatile nature allows it to coordinate the recruitment, activation and exchange of many proteins whose combined activities allow for the protection and propagation of eukaryotic genomes (Maréchal and Zou, 2014). How multiple RPAs associate on ssDNA and coordinate its vast array of processes remains to be determined (Sun et al., 2019). However, a critical feature of RPA is that, though it can bind nucleic acids with very high affinity, it can easily be displaced by other enzymes for further downstream processing (Sun et al., 2019).

When ssDNA length is sufficient for HR repair, the end resection process is terminated (Zhao et al., 2015). Although the regulation of DNA end resection termination are not yet clearly understood, several studies suggest that under physiological conditions, end resection is terminated by RAD51-RPA switching (Zhao et al., 2015). This switching is regulated by BRCA2-DSS1. DSS1 - SEM1 in yeast - is a small, highly acidic protein that competes with ssDNA, by mimicking ssDNA in order to remove RPA from the genuine ssDNA (Zhao et al., 2015; Stefanovie et al., 2019; Le et al., 2020; Rossi et al., 2021). The DSS1 then binds to BRCA2 in order to facilitate RAD51 filament formation (Liu et al., 2010; Stefanovie et al., 2019; Rossi et al., 2021). DSS1 does not seem to bind DNA on its own but appears to enhance ssDNA binding activities of BRCA2 and RAD52 to promote DSB repair (Zhao et al., 2015). BRCA2 then recruits RAD51 to complete the switch (Zhao et al., 2015; Rossi et al., 2021). For cells with DNA already replicated prior to cell division, RAD51 will oligomerise and form a nucleoprotein filament on the resected, single-stranded DNA (ssDNA) end of the DSB, and search for the homologous DNA sequence on the undamaged sister chromatid, performs strand exchange (invasion), and produce a joint molecule called a D-loop (Rossi et al., 2021). DNA polymerase will then use the homologous DNA strand as a template from the D-loop, and the 3′-end of the broken DNA strand as a primer to commence DNA repair synthesis (Rossi et al., 2021). The other end of the double-strand break is then apprehended by RAD52, joining it to the D-loop, through the annealing process, causing the displaced strand to act as a template for the second strand synthesis (Rossi et al., 2021). When DNA synthesis is complete, the D-loops are then dissociated by RAD54, a protein that interacts with RAD51 to promote branch migration, or interacts with helicases like BLM (van Brabant et al., 2000; Bugreev et al., 2006; Kawale and Sung, 2020; Rossi et al., 2021). DNA is further extended by DNA polymerase, annealed to the ssDNA part of the second broken DNA, gap filled and finally restored (Rossi et al., 2021). In Figure 1 we summarise what is described in this section about the DSB signalling mechanisms.

Homologous recombination is able to repair DSBs error-free using the undamaged sister chromatid (Dietlein and Reinhardt, 2014). As the accuracy of homologous recombination repair is important for DSBs (Sugiyama and Kantake, 2009), if it is impaired, chemotherapeutic opportunities may arise (Huang and Zhou, 2021).

We considered HDR through Homologous Recombination (HRR) network as available in Pathways Commons (Cerami et al., 2010) in the SIF (Simple Interaction Format) format at the link in the reference (Orlic-Milacic, 2015). See these data also reported in Supplementary Tables S1–S3.

We implemented an R script, that takes as input the HR network and is able to.

• analyse the topology of the network through the calculation of standard and new centrality measures. The standard node centrality measures considered in this study are the degree (in-, out- and total), the betweenness, the clustering coefficient, the eingenvector centrality, the vibrational centrality, the subgraph centrality, and the information centrality (see (Marsden, 2005; Koschützki and Schreiber, 2008; Ghasemi et al., 2014; Wang et al., 2014; Fornito et al., 2016; Jalili et al., 2016; Ashtiani et al., 2018) for a concise but comprehensive report on the meaning and the use of these measure in molecular biology). We considered also a new centrality measures, such as vibrational centrality, introduced by Estrada in (Estrada and Hatano, 2010) that we will discuss in more detail in the next section (we also refer the reader to (Lecca and Re, 2019) for a review on vibrational centrality); for the reader’s convenience, we list the definition of these centrality indices in Supplementary Table S4, that are also extensively covered in many textbooks on graph theory, and in various articles in the applied sciences. We refer the reader to Estrada’s numerous works, a comprehensive compendium of which can be found in the book (Estrada and Hatano, 2010; Estrada, 2011);

and a R script implementing parametric sensitivity analysis of the dynamics model.

In the dynamics model, by default the kinetic rate constants k a well as the initial values of the proteins and molecules concentrations are set equal to random values in fixed ranges. Nevertheless these ranges can be modified by the user as shown by the interactive console output reported in Table 1. However, we note that the interval of definition of the uniform distribution cannot exceed the maximal range of parameter variability within which the system of rate equations has a solution. We refer the reader to a previous work of us (Lecca et al., 2016) for more details on this.

In a first experiment, the initial values of the proteins concentration (not experimentally known) has been drawn randomly in a range [1, 100] a. u., and the simulation time interval was [0, 10] a. u. In a second experiment, the numerical simulation of the model was performed by assigning an initial quantity between 18 and 20 (expressed in arbitrary units a. u.) to each node and for t ∈ [0, 1400] (in arbitrary units). The solution of the system of 25 differential equations converges for values of rate constants in the range [ ] 0,0.01 a. u. in the first experiment, and in the range [ ] 1 0 − 6 , 1 0 − 5 a.u. in the second experiment. To the best of our knowledge, the in vivo concentrations of the proteins that are part of the HR network are not known. Combined with the lack of knowledge of the values of rate constants, this means that we cannot assign units to these values that reflect in vivo kinetics. However, the initial values of the protein concentrations have been chosen from the order of magnitude at which various in vitro experiments operate and express these concentrations in nM and time in seconds. The literature we have referred to includes the works of (Yang et al., 2010; Nguyen et al., 2013; Foertsch et al., 2019).

We changed the parameters one at a time while keeping the values of the others fixed. Since for each parameter p _h (h ranges from 1 to the number of parameters in the equations), we sampled N _P values, and consequently we performed N _P model simulations. Let us denote with x _s (t), (s = 1, 2, …, d) the time series expressing the numerical solutions of the rate equations, where d is the number of the proteins in the network. The index of sensitivity of x _s (t) with respect to the change of h-th parameter from the value p _h to the value p h ′ is calculated as in (Lecca et al., 2016) by S I i s h = 1 N ∑ k = 1 N 1 N P − 1 ∑ r = 1 N P x s r t k | p h ← p h ′ − x s ̄ r p h ← p h ′ 2 1 2 (1) where N is the length of the time series x s (t) , and x s ̄ p h ← p h ′ = 1 N P ∑ k = 1 N P x s t k | p h ← p h ′ . (2) where “ p h ← p h ′ ” means “p _h replaced by p h ′ ”. The Eq. 1 defines the mean of the standard deviations of the distributions of the simulated values of a protein/gene abundance at time points t _k .

In case the user knows the values of the rate constants, he/she can add them as an extra column to the SIF format of the input files. At the moment of writing, for most interactions the values of the kinetic constants are not known and there are no time-resolved data from which it is possible to infer them. It is precisely this context that justifies our choice to study the dynamics of the system in a range of values of the model parameters and more generally to provide a software that can be used as a platform for in silico experiments.

We implemented Network Analyser and Dynamics Simulator (NADS) consists of three modules written in R language:

• network_analysis_functions.R: this module implements the functions that processes the SIF data-frames to make them suitable to their conversion into a graph. This module also implement the functions to rank the nodes according to their centrality measures;

The user can launch the software simply by running in RStudio the script network_analysis.R, and then by answering to the questions in the interactive interface as shown in Table 2.

We provide also the fourth module implementing the parametric sensitivity analysis, named sensitivity_analysis.R, which takes as an input the system of equations automatically generated by network_analysis.R.

The main module is the script network_analysis.R. As soon as the user runs it from the R Studio GUI (R Studio, 2022) or from a terminal, an interactive output is displayed as in Tables 2, 3. The program asks the user to select the network to be analysed and then.

• it calculates the centralities measures

The program returns also the execution times for the tasks expected to be the most computationally demanding, such as integrating the equation.

The HR network considered in this study has 25 nodes and 250 edges, as reported in Supplementary Tables S1–S3. The right part of the Figure 4 shows the HR network in circular layout. The network analysis phase of our study calculated the centrality measures distributions show in Figure 3, and identified six genes, as shown in Table 4:

1. BLM, scoring first for vibrational centrality

Of particular interest is the fact that BLM has the highest vibration centre. The interpretation of this result is that BLM is the node most sensitive to stresses and/or stimuli, i.e., according to the vibrational centrality measure, it is the most vulnerable node in the network (Estrada and Hatano, 2010). This result is of particular interest in light of the crucial role this gene plays in the HR network. Indeed, the key role of BLM is well know, and alterations in this protein is linked to different diseases including cancer (Kaur et al., 2021). BLM is a 3′-5′ ATP-dependent RecQ DNA helicase. It is a genome stabilizer playing an essential role in the DNA replication regulation, DNA recombination, and both homologous and non-homologous pathways of DSB repair. The high vulnerability of the BLM node to external stimuli and conditions suggests the need to identify which conditions and/or stimuli may be altering it, in order to preserve its proper functioning and/or to understand how it can be restored if it is altered. The high vulnerability of this node could also be explained by a recent study by Kaur et al. (2021). These authors report that BLM has a dual function both as a tumour suppressor and possibly as a proto-oncogene, being probably involved in the mechanisms of its deregulation in tumours.

The analysis also correctly identifies the SEM1 gene as a node with a high clustering coefficient. Indeed, as reported in (Safran et al., 2021; GeneCard, 2022), SEM1 gene encodes for a protein that is part of a 26S proteasome, which is a multiprotein complex with a function in the ATP-dependent degradation of ubiquitinated proteins. This complex contributes to the maintenance of protein homeostasis by removing misfolded or damaged proteins, which could jeopardize the healthy cellular functions, and by removing proteins no longer need. Therefore, 26S proteasome is involved in numerous cellular processes, including cell cycle progression, apoptosis, or DNA damage repair (Sone et al., 2004).

SEM1 was found also as a subunit in experiments of affinity purification of the yeast 19S proteasome, and its human homolog, DSS1, was found to copurify with the human 19S proteasome (Krogan et al., 2004). DSS1 is associated with the tumour suppressor protein BRCA2 involved in DNA DSBs repair. The authors in (Krogan et al., 2004) proved that SEM1 is essential for efficient repair of an HO-generated yeast DSB using both HR and nonhomologous end joining (NHEJ) pathways. Moreover, they showed that deletion of SEM1 contributes to cause defects in (synthetic) growth and hypersensitivity to genotoxins when combined with mutations in certain well-established genes involved in the DNA DSB repair.

Similarly to SEM1, the result of a high clustering coefficient is also expected for RPA4, as RPA4 is also part of a complex (Keshav et al., 1995). RPA4 gene encodes a single-stranded DNA-binding protein that is a subunit of the replication protein A complex (GeneCard, 2022). Replication protein A is essential for DNA DSB repair and plays a crucial role in the activation of cell cycle checkpoints. As regards the RPA complex, we have already seen in the previous sections that the RPA complex controls DNA repair and DNA damage checkpoint activation as well. In particular, the network analysis shows that RPA1 highly scores by total degree and betweenness. These results reflect the fact that RPA1 is an active route of communication exchanges between various nodes in the network. RPA1 is part of the heterotrimeric replication protein A complex (RPA/RP-A). It stabilizes single-stranded DNA intermediates, that form during DNA replication or upon DNA stress In (Bass et al., 2016; Haahr et al., 2016; Human Protein Atlas, 2022). It prevents the reannealing of single-stranded DNA intermediates and recruits and activates different proteins and complexes forming part of DNA metabolism. Thereby, it is a key protein both in DNA replication and in the cellular response to DNA damage (Lin et al., 1998). RPA2 shows high score for in-degree and eigenvevtor centrality, meaning that it is interacting with protein also highly scoring by eigenvector centrality and degree (Hansen et al., 2020), and thence with proteins which have a great influence in the HR network. Indeed RPA 2 gene has been found highly expressed in low grade carcinomas and its expression has a gradual significant decrease from stage I to stage IV carcinomas. All the three subunits RPA1, RPA2, and RPA3, were more abundant (with statistical significance evidence) in lymph node negative and earlier stage (stage I and II) gastric carcinomas (Fourtziala et al., 2020). Finally, of particular interest and the fact that RPA3 ranks first in terms of centrality out-degree, hub-centrality and sub-graph centrality. Since subgraph centrality of a node is the number of subgraphs a node participates in (weighted according to their size) (Estrada and Rodríguez-Velázquez, 2005), it means that RPA3 take part into a number of subgraphs of significant size relatively to the whole network size. From our analysis it results that RPA3 versus RPA1 and RPA2, although its roles are similar to those of RPA2, has thus a great influence on pathways of crucial importance more than on subset of unconnected nodes or single nodes.

The analysis also highlights RAD50 that is a component of the MRN complex The protein complex is involved in numerous enzymatic activities required for nonhomologous joining of DNA ends. It is protein is essential for DNA double-strand break repair, cell cycle checkpoint activation, telomere maintenance, and meiotic recombination. (Carney et al., 1998; de Jager et al., 2001; Estrada and Ross, 2018; Bian et al., 2019; Beikzadeh and Latham, 2021; National library of Medicine, 2022). This role is reflected by th ehigh sungraph centrality that measures the centrality of a node by taking into account the number of subgraphs the node participates in. Specifically, the subgraph centrality of a node is the number of closed loops originating at the node, where longer loops are exponentially downweighted. Consequently, subgraph centrality measures how close a node is to the other nodes in the network.

The results of this analysis reveal a correspondence between the measure of centrality and the role of the protein. On the basis of this, when the role of the protein is known, this information can be used to work out the correctness of the computational analysis. When, on the other hand, the role of the protein is not known, knowledge of its centrality measurements can suggest the type or set of types of possible roles.

In Tables 5, 6 we show the rate equations of the HR network dynamics generated as the automatic translation of the network (see the script dynamics.R in the GitLab repository of NADS software). In file Simulations_of_Dynamics_HR_pathway.pdf provided in the Supplementary Material, we show the time evolution curves of each node of the HR network obtained as a solution of the equations.

The parameter sensitivity analysis was conducted by perturbing each parameter in the convergence range of the solution and yielded the results shown in Table 7; Figures 5, 6. In Table 7 we report the kinetic rate constants for which the sensitivity index belongs ot the 98th percentile of the sensitivity index distribution. They correspond to the most sensitive parameters, i.e. to the interactions whose alterations can significantly alter the dynamics of the network. To find out which interactions they refer to, the reader can refer to the Supplementary Tables S1–S3. Figure 5 we shows that MRE11, followed by POLD1, BLM has the highest average sensitivity index. In Figure 6 we show the coefficient of variation of each protein in the HR network. The coefficient of variation, being the ratio of the standard deviation to the mean, measures the extent of variability in relation to the mean of the population. The higher the coefficient of variation, the greater the dispersion. BLM, followed by MRE11 and RPA2 exhibits the lowest coefficient of variation of the sensitivity index. This mean that BLM, MRE11 and RPA2 have high sensitivity indices, and that the distribution of the sensitivity indices is well shaped around its mean, i.e. these protein exhibit almost the same sensitivity for all the parameters of the model. In this study we have therefore found that the BLM and RPA2 are sensitive nodes and that their sensitivity has two components: a topological sensitivity expressed by vibrational centrality, eigenvector centrality and clustering coefficient, and a dynamic sensitivity expressed by the parameter sensitivity index.

As with the BLM and RPA2 proteins, sensitivity analysis also highlights the MRE11 protein, which is highly sensitive to kinetic parameters, and its vibrational centrality is 1.0587524, very close to the maximum value exhibited by BLM (see Table 4). Its eigenvector centrality is 0.7397 which, although not the maximum, is very close to it (see Table 4). Indeed, MRE11 is an integral part of the protein complex of RAD50-MRE11A-NBS1 known as the MRN complex (Porras, 2014; Shibata et al., 2014; Mukherjee et al., 2019). It plays a key role in homologous recombination, and it is generally believed that MRE11 initiates double-strand breaks resection. In particular, the authors show that the loss of MRE11 reduces the efficiency of homologous recombination in human TK6 cells without affecting double-strand breaks resection, indicating a role for MRE11 in homologous recombination also at a post-resection step.

The high value of the eigenvector centrality fork BLM, RPA2 and MRE11 confirms the crucial role of these proteins in the network and expresses the fact that they are pointed by nodes that have a high value of the eigenvector centrality too. Indeed, if a node is pointed to by many nodes (which also have high eigenvector centrality) then that node will have high eigenvector centrality (Fletcher and Wennekers, 2018). The high sensitivity to the parameters characterising the dynamics of the interactions between these proteins and the partners pointing to them indicates the great influence that these partner nodes have on these proteins. Interestingly, RAD51 does not result sensitive to any parameter. The RAD51 encodes a protein that is essential for repairing damaged DNA. Recent findings have indicated RAD51 protein is overexpressed in a variety of tumours Chen et al. (2017). The overexpression of RAD51 causes improper and hyper-recombination, and thus contributes to genomic instability and genetic diversity. Genomic instability might, in turn, drive regular cells towards neoplastic transformation or further contributes to cancer metastatic progression (Chen et al., 2007). The RAD51 protein binds to the DNA at the site of a break and encapsulates it in a protein sheath, initiating the repair process MedlinePlus (2022); Uniprot (2022). RAD51 protein interacts with BRCA1 and BRCA2, to fix damaged DNA. The BRCA2 protein regulates the activity of the RAD51 protein by transporting it in the nucleus to sites of DNA damage. Although the interaction between the BRCA1 protein and the RAD51 protein has still to be elucidate, research suggests that BRCA1 may also activate RAD51 in response to DNA damage (Cousineau et al., 2005; Chappell et al., 2016). The result of the sensitivity analysis found seems to contradict the important role of this protein in these interactions. Indeed, for example, one might expect a high sensitivity of RAD51 to the k139 due to its interaction with BRCA2 (see Supplementary Table S2). One explanation for this contradiction could be that since the mechanisms of interaction of RAD51 with these proteins are not fully known, the model used in this study could be an oversimplification of the interaction of RAD51 with its partners. If more accurate models in the future confirm the low sensitivity of RAD51 to the parameters of the rate equations describing the dynamics of the network, it will be necessary to investigate the physical and biological characteristics that make it so stable to perturbations. The fact that RAD51 has a low value of vibrational centrality in this study is a factor in favour of the possible confirmation of this case.

In Figures 7, 8 we report the results of the sensitivity analysis obtained selecting different ranges of initial conditions and parameters values. The plots highlights RAD51C, MRE11, RAD50 and BRCA2 as the most sensitive nodes to the parameters. We comment in the Section Remarks the expected differences and similarities in the results of sensitivity analysis when we change the intervals of initial conditions and parameters.

This report presented an application of network analysis and mathematical modelling to the double-strand break repair pathway homologous recombination repair (HR). The complexity of the network of repair mechanisms itself, as well as the complexity of its interactions with the surrounding environment (Li et al., 2009; Chatterjee and Walker, 2017; Kusakabe et al., 2019; Poetsch, 2020; Roux et al., 2021), and the mutations of its components make its mathematical modelling particularly difficult, especially when based on rate equations. It is therefore of great necessity to have a tool that can implement these two important steps:

1. network analysis including standard centrality measures and new measures to quantify the robustness and responsiveness of the network to stimuli and stresses not dependent on the network topology

The implementation of these step is an innovative perspective for the analysis of DNA repair mechanisms. So far in the literature, there are many studies and analyses focused on the genetic and genomic aspects of the pathways, but studies on the mathematical modelling of its dynamics are absent. Our study therefore aims to fill this gap, since the mechanisms of DNA repair are governed by genes, proteins and pathways in continuous communication with the environment. For this reason, the analysis of the dynamics of the network is particularly useful, since it can quantify the vulnerability of the network and the modes of response to stimuli and exogenous stress. To the best of our knowledge there are no schemes for translating a graph associated with a biological network into a set of dynamic equations. The main reason for this is that there is no unambiguously defined semantics of a graphic representation of a biological network, i.e. there is no unambiguous definition of the graphic symbolism in terms of the mathematical equation describing the interaction indicated by that symbolism. The translation model we propose in this study is a basic model that describes the interactions indicated by the graph with linear first-order differential equations in the parameters. The code that implements this translation, however, gives the user the possibility to modify the model where he/she deems it appropriate in the light of available biological knowledge, or in the case he/she like to generate new hypothetical scenarios. We believe that the availability of a tool such as NADS can support the investigation of such a complex network that is subject to continuous interaction with external agents, not only to understand its dynamics, but also to predict its evolution and identify points of vulnerability for the benefit of the medical applications that this research may provide.