Inequality in the Frequency of the Open States Occurrence Depends on Single 2H/1H Replacement in DNA

In the present study, the effect of 2H/1H isotopic exchange in hydrogen bonds between nitrogenous base pairs on occurrence and open states zones dynamics is investigated. These processes are studied using mathematical modeling, taking into account the number of open states between base pairs. The calculations of the probability of occurrence of open states in different parts of the gene were done depending on the localization of the deuterium atom. The mathematical modeling study demonstrated significant inequality (dependent on single 2H/1H replacement in DNA) among three parts of the gene similar in length of the frequency of occurrence of the open states. In this paper, the new convenient approach of the analysis of the abnormal frequency of open states in different parts of the gene encoding interferon alpha 17 was presented, which took into account both rising and decreasing of them that allowed to make a prediction of the functional instability of the specific DNA regions. One advantage of the new algorithm is diminishing the number of both false positive and false negative results in data filtered by this approach compared to the pure fractile methods, such as deciles or quartiles.


Introduction
The deuterium concentration in the body plays an important role in the metabolic processes of living systems [1][2][3][4]. The biological effects caused by deuterium depleted water (DDW) have been studied at various levels: cellular [5][6][7][8], tissue [9], organismic [10][11][12][13][14]. It was found that low concentrations of deuterium in the drinking diet increase stress resistance in mammals [15], and may also have neuroprotective properties [16]. In addition, it was shown that a decrease in the deuterium content in the blood and brain of laboratory animals reduces the disruption of antioxidant enzymes functioning during hypoxia in comparison with animals that received water with a natural deuterium content [17]. The administration of water with a reduced deuterium content into the drinking diet leads to a change in the 2 H/ 1 H ratio in the body tissues due to the hydrogen isotopes exchange [18], in particular, as shown by the recent study, the antiproliferation effects of DDW are associated with the occurrence of an imbalance in the mitochondria between the production of reactive oxygen species and their neutralization [19]. Some studies have shown an increase in H 2 O 2 production by mitochondria during their pre-incubation in the deuterium depleted media [20,21]. In addition, the mechanisms of DDW influence on biological systems from the point of view of the medical biochemistry are logically explained by a modification in the mitochondria energy balance [22,23].
Thus, to date, various biological effects caused by DDW have been examined in sufficient detail, and there is also an opinion that the main "fast" mechanism of the effect of low concentrations of deuterium in drinking water on the organism is a change in the parameters of the mitochondria work. However, it is important to understand the long-term effects of DDW administration, as well as its possible impact on DNA [24]. It is known that the deuterium bond energy is~134 calorie /mol higher than the protium energy; therefore, deuterium forms a 5% stronger hydrogen bond [25]. The infiltration of a deuterium atom instead of protium into the hydrogen bond of the DNA double helix can cause a very noticeable temporary malfunction in the transmission of this information, most likely due to a delay in the opening of any hydrogen bond. Deuterium atoms get into the hydrogen bonds of DNA double helices because of the rapid protium-deuterium isotope exchange with the surrounding water molecules [26]. With a natural deuterium content of 156 ppm in water, the equilibrium probability of deuterium getting into each of the possible hydrogen bonds is small and equals approximately 2 × 10 −4 [24]. In addition to fluctuation twists of the helix axis and rotations of adjacent pairs of nitrogen bases in the DNA molecule, the opening and closing of individual pairs of nitrogen bases occurs [27]. This process leads to changes in conformation and plays an important role in the reactions of DNA with chemical agents. Replacing a protium atom with deuterium affects the process of opening base pairs by increasing the energy required to break the bond. To date, the implementation of an experiment that would allow us to evaluate the effect of a deuterium atom on the processes of opening base pairs is difficult. A theoretical study of these processes is possible using methods of mathematical modeling, and one of the key conditions for the adequacy of the mathematical model of DNA is to take into account open states. Existing mathematical models of DNA are discussed in a number of reviews [27][28][29]. However, to describe the effects caused by the introduction of deuterium atoms to hydrogen bonding between base pairs, the upgraded Yakushevich model was chosen [30][31][32][33], which allows to take into consideration the energy of the hydrogen bonds between pairs of nitrogenous bases. Due to the introduction of an additional term, this model allows to consider the effects of dissipation caused by the viscosity of the medium surrounding the DNA molecule. In addition, in the model, due to the introduction of the torque constant of a given site of the sugar-phosphate chain, the interactions between adjacent base pairs are indirectly taken into account. Another argument in favor of choosing the Yakushevich mechanical model to describe the role of deuterium atoms in the opening and closing processes of base pairs was the fact that radial torsion models are able to take into account distortions of the DNA structure caused by external torsion stress. However, most physicochemical experiments are carried out on relaxed DNA, without such distortions. We can assume that the behavior of relaxed molecules is described equally well by radial and radial-torsion models [27].
In the presented work, in the framework of the mechanical DNA model, evaluation experiments were carried out on the effect of isotopic 2 H/ 1 H exchange in hydrogen bonds between pairs of nitrogen bases on the uneven probability distribution of open states along the gene length.

Mathematical Model
To simulate the processes of formation and dynamics of open states (OS) we will use a mathematical model that describes the rotational movement of nitrogen bases around the sugar-phosphate chain of a DNA molecule [34].
(2"') here: ϕ i j (t)-is the angular deflection of the i-th nitrogen base of the j-th chain counted counterclockwise at time t; I i j -is the rotational inertia of the i-th nitrogen base of the j-th chain; R i j -is the distance between the center of inertia of the i-th nitrogen base of the j-th chain to sugar phosphate chain; K i j -is the constant characterizing the torsion moment of the i-th segment of the j-th sugar phosphate chain; k i 12 -is the constant characterizing the bond elastic properties of the i-th nitrogen base pairs; F i j (t)-external influence on the i-th nitrogen base of the j-th chain at a time t, n-is the number of nitrogene base pairs in the system.
The magnitude of the external force is taken as F i j dt (t) models the effects of dissipation due to the interaction with water surrounding the DNA molecule, the term F 0 cos ωt is the external periodic impact.
In Equations (1 )-(2"'), the first term to the right of the sign of the equality describes the force acting on the i-th nitrogen base from the sugar-phosphate chain, the second term is the force from the complementary nitrogen base, and the third term is the external impact.
Thus, Equations (1 )-(2"') allow us to simulate the hydrogen bond in the i-th pair (δ i = 1, k i 12 = k H,i 12 ), deuterium bond (δ i = 1, k i 12 = k D,i 12 ) and disruption of this bond (δ i = 0). We will assume that the break between base pairs occurs if the potential binding energy in these pairs exceeds a certain critical value E H cr for hydrogen bonds and E D cr for deuterium bonds, but if the potential energy in a pair with a broken bond is less than the critical value, then the bond is restored.
We add the initial conditions to Equations (1 )-(2"'): For definiteness, we'll assume that at t = 0 the system is in equilibrium, i.e., in the initial conditions (3) Problem (1)-(3) is the Cauchy problem for a system of 2n ordinary differential equations; in this work, all studies were carried out on the basis of a numerical solution of this system.
Problem (1)-(3) is the Cauchy problem for a system of 2n ordinary differential equations; in this work, all studies were carried out on the basis of a numerical solution of this system.
On Figure 1

The Effect of 2 H/ 1 H Exchange on the Probability of OS Formation
We will study the effect of 2 H/ 1 H exchange on the formation and dynamics of the OS using the example of a gene encoding interferon alpha 17. For this gene n = 980, the values of the Equations (1′)-(2‴) coefficients are shown in Table 1 (data taken from [30]), F0 = 0.526 • 10 −22 J, ɷ = 0.4 • 10 12 s −1 . We designate by P0 the probability of an OS formation in a molecule in which all pairs of nitrogenous bases are connected by hydrogen bonds; by Pi, I = 1, , the probability of an OS occurrence in a DNA molecule in which in the i-th nitrogen base pair one <any > hydrogen bond is replaced by deuterium.
The probabilities P0 и Pi, I = 1, , will be sought on the basis of a numerical solution of the problem (1)-(3). To do this, we will create a set of points tj = jƮ, j = 1, m, Ʈ = T/m in the segment [0, T]. Then, we will calculate at t = tj the ratio qj of the number of base pairs with a broken bond to the total

The Effect of 2 H/ 1 H Exchange on the Probability of OS Formation
We will study the effect of 2 H/ 1 H exchange on the formation and dynamics of the OS using the example of a gene encoding interferon alpha 17. For this gene n = 980, the values of the Equations (1 )-(2"') coefficients are shown in Table 1 (data taken from [30]), F 0 = 0.526 × 10 −22 J, ω = 0.4 × 10 12 s −1 . We designate by P 0 the probability of an OS formation in a molecule in which all pairs of nitrogenous bases are connected by hydrogen bonds; by P i , I = 1, n, the probability of an OS occurrence in a DNA molecule in which in the i-th nitrogen base pair one <any > hydrogen bond is replaced by deuterium.
The probabilities P 0 and P i , I = 1, n, will be sought on the basis of a numerical solution of the problem (1)- (3). To do this, we will create a set of points t j = jτ, j = 1, m, τ = T/m in the segment [0, T]. Then, we will calculate at t = t j the ratio q j of the number of base pairs with a broken bond to the total number of base pairs n, then the value of P k is equal to the arithmetic mean value over the points t j of these ratio: Since the deuterium bond is 5% stronger than hydrogen [25], we took For various values of E H cr values P 0 and P i , i = 1, n were calculated. The calculations were performed with accuracy 10 −6 for T = 3.0 × 10 −10 s, τ = 0.0001 × 10 −10 s.
To carry out calculations based on the mathematical model, a computer program was developed. It works as follows: when a deuterium atom enters hydrogen bonds between the first base pair of gene encoding interferon alpha 17,980 equations are solved (by the Runge-Kutta method of the 4th order), the average probability of occurrence of open states in the gene is calculated, then a similar calculation is performed for the second pair of bases, etc.

Results
The OS dynamics of the DNA molecule and the effect of 2 H/ 1 H exchange on it are illustrated on

Results
The OS dynamics of the DNA molecule and the effect of 2 H/ 1 H exchange on it are illustrated on      Table 2) and black dots represent coincidence values for i min and i max . Table 2 contains data of the values of P 0 , the i min and i max numbers of nitrogen base pairs, the substitution of which in the hydrogen bond with deuterium bond leads to the lowest and highest values of the open state occurrence probability and the values of the P i min and P i max . Gene encoding interferon alpha 17 (IFNA17) and containing 980 nucleotide pairs were conditionally divided into three equal parts: I part (from the 1st to the 327th bases, n = 327), II part (from the 328th to the 653th bases, n = 326) and III part (from the 654th to the 980th bases, n = 327). In the whole gene we counted probabilities (P i ) of occurrence of open states between different nitrogenous bases in double-stranded DNA dependent on the single 2 H/ 1 H replacement in base pair of each gene region. All these P i were arranged from P imin to P imax and each of them was compared to P 0 , which was determined as probability of OSs occurrence, when all hydrogen bonds in DNA are 1 H. According to the serial number of each nucleotide base pair its affiliation to I part, II part or III part of IFNA17 was determined, and after that the amount of nucleobase pairs from the ranges "Maximum" and "Minimum" for the whole gene and its parts was calculated. Nitrogenous bases, which had P imax , and base pairs with the higher P i were selected as representatives of the maximum range ("Maximum"). Nitrogenous bases, which had P imin , and base pairs with the more-low P i were selected as representatives of the minimum range ("Minimum"). The ranges "Maximum" and "Minimum" were found by different below-presented methods: decile method [35], quartile method [35], and new method (Basov-Jimack algorithm, or BJ-algorithm).
The range of E H cr was from 0.30 to 0.59, while the change in energy taken into account when calculating open states in the entire range was uniform and equaled to 0.01 (Table 2).
In addition, to exact the number of n min and n max of the nitrogen base pairs in the range from 0.30 to 0.59 (where each E H thr had n min or n max more than 0 for the whole gene), their number was also calculated for smaller ranges of E H cr : a.
from 0.30 to 0.43, which were selected because of III part, which, in this range, had summary n max more than 0; its n max measured by decile method equaled 203 (Figure 3), quartile method equaled 88 (Figure 4), or BJ-algorithm equaled 87 ( Figure 5); b. from 0.44 to 0.59 that were selected because of III part, which in this range had summary n max equaled 0; its n max , measured by the each approach: decile method, quartile method, or BJ-algorithm, always equaled 0 ( Figure 3). a. from 0.30 to 0.43, which were selected because of III part, which, in this range, had summary nmax more than 0; its nmax measured by decile method equaled 203 (Figure 3), quartile method equaled 88 (Figure 4), or BJ-algorithm equaled 87 ( Figure 5); b. from 0.44 to 0.59 that were selected because of III part, which in this range had summary nmax equaled 0; its nmax, measured by the each approach: decile method, quartile method, or BJ-algorithm, always equaled 0 ( Figure 3). After that, we used some statistical methods to exact the significance of differences among nmin and nmax of the three gene parts. Yates corrected chi-squared test (χ 2 Yates) was applied for a 2 × 2 contingency Table 3 (where degrees of freedom (ν) = (2-1)•(2-1) = 1). As a shortcut, for a 2 × 2 Table 3 with the following entries: Table 3. Chi-squared test with Yates correction.

S F
Where A and B are the rows according to the parts of the gene: I, II, or III; S is the column of nmin of the range "Minimum", or nmax of the range "Maximum"; F is the column of the rest number of nucleotide pairs; a and c are nmin of the range "Minimum", or nmax of the range "Maximum" in each compared part of the gene in the determined range of E H thr; b and d equals the total number of nucleotide pairs of the gene part minus nmin of the range "Minimum", or nmax of the range After that, we used some statistical methods to exact the significance of differences among n min and n max of the three gene parts. Yates corrected chi-squared test (χ 2 Yates ) was applied for a 2 × 2 contingency Table 3 (where degrees of freedom (ν) = (2-1)·(2-1) = 1). As a shortcut, for a 2 × 2 Table 3 with the following entries: Table 3. Chi-squared test with Yates correction.

S F
Where A and B are the rows according to the parts of the gene: I, II, or III; S is the column of n min of the range "Minimum", or n max of the range "Maximum"; F is the column of the rest number of nucleotide pairs; a and c are n min of the range "Minimum", or n max of the range "Maximum" in each compared part of the gene in the determined range of E H thr ; b and d equals the total number of nucleotide pairs of the gene part minus n min of the range "Minimum", or n max of the range "Maximum" in each compared part of the gene in the determined range of E H thr ; "Maximum" in each compared part of the gene in the determined range of E H thr; NA = a + b; NB = c + d; NS = a + c; NF = b + d; N = NA + NB + NS + NF. Chi-square corrected by procedure for Bonferroni (χ 2 B) was used for a 3 × 2 contingency table (where 3 is rows, 2 is columns, ν = (3-1)•(2-1) = 2). The Kruskal-Wallis ANOVA by Ranks test (KWt) was applied for the comparison nmin of the range "Minimum", or nmax of the range "Maximum" in the determined range of E H cr for two gene parts that are mutually independent. Median test (Mt) was used for the comparison nmin of the range "Minimum", or nmax of the range "Maximum" in the determined range of E H cr for two gene parts that are mutually independent.
It was found out that the quantity of nmin and nmax of OSs for each gene part was different in the above-described ranges of E H cr (Table 4). In the range of E H cr from 0.30 to 0.59 the highest number of nmax was in the II part, e.g., when it was counted by the decile method (χ 2 B: p < 0.001) the nmax in the II part was in 2.76 and 2.89 times higher than in the I part and III part respectively, but for the quartile method (χ 2 B: p < 0.001) the nmax in the II part was in 5.35 and 6.69 times higher, and according to BJalgorithm (χ 2 B: p < 0.001) the nmax in the II part was in 4.74 and 5.94 times higher compared to the I part and III part respectively (Table 4). Wherein in the same range of E H cr in the II part of IFNA17, the highest number of nmin also was only if it was counted by BJ-algorithm (χ 2 B: p < 0.002), which showed that the nmin in the II part was in 1.22 and 1.30 times higher than in the I part and III part, respectively. Whereas the quantity of nmin calculated by the quartile method turned higher in the I part compared to the II part and III part by 9.2% and 48.5% respectively (χ 2 B: p < 0.001, Table 4). It should be noted, that nmin, counted by the decile method at the same condition, had no differences in all parts of the gene (χ 2 B: p > 0.2). Chi-square corrected by procedure for Bonferroni (χ 2 B ) was used for a 3 × 2 contingency table (where 3 is rows, 2 is columns, ν = (3-1)·(2-1) = 2). The Kruskal-Wallis ANOVA by Ranks test (KWt) was applied for the comparison n min of the range "Minimum", or n max of the range "Maximum" in the determined range of E H cr for two gene parts that are mutually independent. Median test (Mt) was used for the comparison n min of the range "Minimum", or n max of the range "Maximum" in the determined range of E H cr for two gene parts that are mutually independent.
It was found out that the quantity of n min and n max of OSs for each gene part was different in the above-described ranges of E H cr (Table 4). In the range of E H cr from 0.30 to 0.59 the highest number of n max was in the II part, e.g., when it was counted by the decile method (χ 2 B : p < 0.001) the n max in the II part was in 2.76 and 2.89 times higher than in the I part and III part respectively, but for the quartile method (χ 2 B : p < 0.001) the n max in the II part was in 5.35 and 6.69 times higher, and according to BJ-algorithm (χ 2 B : p < 0.001) the n max in the II part was in 4.74 and 5.94 times higher compared to the I part and III part respectively (Table 4). Wherein in the same range of E H cr in the II part of IFNA17, the highest number of n min also was only if it was counted by BJ-algorithm (χ 2 B : p < 0.002), which showed that the n min in the II part was in 1.22 and 1.30 times higher than in the I part and III part, respectively. Whereas the quantity of n min calculated by the quartile method turned higher in the I part compared to the II part and III part by 9.2% and 48.5% respectively (χ 2 B : p < 0.001, Table 4). It should be noted, that n min , counted by the decile method at the same condition, had no differences in all parts of the gene (χ 2 B : p > 0.2).  Note: * is p-value (χ 2 Yates) < 0.05 compared to the similar n of the I part, measured by the same method; # is p-value (χ 2 Yates) < 0.01 compared to the similar n of the I part, measured by the same method; ¤ is p-value (χ 2 Yates) < 0.01 compared to the similar n of the II part, measured by the same method.  Note: * is p-value (χ 2 Yates ) < 0.05 compared to the similar n of the I part, measured by the same method; # is p-value (χ 2 Yates ) < 0.01 compared to the similar n of the I part, measured by the same method; ¤ is p-value (χ 2 Yates ) < 0.01 compared to the similar n of the II part, measured by the same method.
Moreover, the comparison of the total n max number of the gene ends (sum of I part and III part) to the base pair amount of the range "Maximum" in the center gene part revealed that the II part had more n max than the both ends of IFNA17; which was 1.41 times more by the decile method (pχ 2 Yates < 0.0001), 2.97 times more by the quartile method (pχ 2 Yates < 0.0001), and 2.64 times more by BJ-algorithm (pχ 2 Yates < 0.0001). In contrast, the nucleobase quantities of the range "Minimum" had no difference between ends and center part of gene calculated by decile method (pχ 2 Yates = 0.057), but n min difference counted by quartile method was significant (pχ 2 Yates = 0.012), and higher than the difference by BJ-algorithm (pχ 2 Yates = 0.0001). The comparing of nucleobase quantities of the range "Maximum" depending on each E H cr (from 0.30 to 0.59) in the different parts of gene showed, that there was the difference among them, and if it was counted by the decile method, by the quartile method, or by BJ-algorithm then the p KWt and p Mt always were less than 0.0001 (Figures 3-5). In contrast of it the nucleobase quantities of the range "Minimum" had no difference when counting by the decile method (p KWt = 0.223; p Mt = 0.325, Figure 3), but the difference was found for n min computed by the quartile method (p KWt = 0.046; p Mt = 0.079, Figure 4) and by BJ-algorithm (p KWt = 0.034; p Mt = 0.111, Figure 5). Additionally, imparity between n max of the ends of IFNA17 and base pair quantity of the range "Maximum" in its center was also proved by the decile method, quartile method, and BJ-algorithm (p KWt < 0.001; p Mt < 0.001 in all cases). However, it did not achieve any significant p KWt or p Mt by all approaches to the n min , and each p KWt and p Mt was more than 0.12 in this range of E H cr .
At the same time in the range of the higher E H cr (from 0.44 to 0.59), the biggest base pair quantity of the range "Maximum" was also in the II part (χ 2 B : p < 0.001), e.g., that n max number calculated by the each method in the II part of the gene was more than 15.5 times higher than in its I part (pχ 2 Yates < 0.0001, Table 4). Moreover, for the III part of IFNA17, there were no nucleotides in the range "Maximum" at all for these energies, so n max of the center part of gene was also more than 15.5 times higher compared to n max of its ends (pχ 2 Yates < 0.0001). Comparison of n max in the three parts of gene, depending on the each E H cr from 0.44 to 0.59, pointed out to the significant difference among them, and if it was counted by the decile method, by the quartile method, or by BJ-algorithm then the p KWt and p Mt always were less than 0.0001 (Figures 3-5). In contradistinction to this, if n min was counted by decile method (χ 2 B : p > 0.5), and by BJ-algorithm (χ 2 B : p > 0.2) the nucleobase quantities in the range "Minimum" had no differences among three parts of IFNA17, but the difference was found by the quartile method (χ 2 B : p < 0.002, Table 4). In addition, there was no statistical inequality of n min between II part and III part (pχ 2 Yates from 0.0816 to 0.1947 dependent on the counting method), ends and center of the gene (pχ 2 Yates from 0.1699 to 0.3313 dependent on the counting method), but if it was computed by the quartile method the base pair amount in the range "Minimum" of the I part was higher compared to the II part and III part by 16.9% (pχ 2 Yates = 0.0099) and 27.1% (pχ 2 Yates = 0.0001) respectively (Table 4). Moreover, when comparing n min in the I, II, and III parts of IFNA17, depending on the each E H cr from 0.44 to 0.59, the significant difference among them was found, which has p KWt = 0.0034 and p Mt = 0.0080 if it was counted by the decile method, p KWt = 0.0006 and p Mt < 0.0001-by the quartile method, and also p KWt = 0.0002 and p Mt < 0.0001-by BJ-algorithm. The comparison of the ends and center of the gene, depending on the each E H cr in this energy range, revealed significant differences between both n max and n min by all calculating methods (for n max : p KWt were less than 0.0001, and p Mt -from less than 0.0001 to 0.0010; for n min : p KWt were from 0.0003 to 0.0060, and p Mt -from 0.0012 to 0.0131).
In the range of E H cr from 0.30 to 0.43 according to the decile method (χ 2 B : p < 0.001, Figure 3) in the II part of IFNA17, the total nucleotide number of the range "Maximum" was higher than n max of the I part in 1.65 times (pχ 2 Yates < 0.0001) and n max of the III part in 1.58 times (pχ 2 Yates < 0.0001). The quartile method demonstrated (χ 2 B : p < 0.001, Figure 4) that in this energy range the II gene part had the total number of n max significant higher compared to n max of the I part (in 3.10 times, pχ 2 Yates < 0.0001), or it of the III part (in 3.24 times, pχ 2 Yates < 0.0001). According to BJ-algorithm (χ 2 B : p < 0.001, Figure 5), the similar number in the II part of IFNA17 was higher than n max of the I part in 2.74 times (pχ 2 Yates < 0.0001) and n max of the III part in 2.89 times (pχ 2 Yates < 0.0001, Table 4). The comparison of the base pair amounts of the range "Minimum" among three parts of the gene showed differences by the decile method (χ 2 B : p < 0.002), quartile method (χ 2 B : p < 0.001), and BJ-algorithm (χ 2 B : p < 0.001). Moreover, gene part, included the highest n min , was different depending on the counting method, e.g., if it was calculated by the decile method the III part had the biggest amount of n min : by 25.3% and 23.9% more than the I part (pχ 2 Yates < 0.0002) and II part (pχ 2 Yates < 0.0005) respectively; and if it was computed by the quartile method then the I part had the highest n min : by 79.4% more than the III part (pχ 2 Yates < 0.0001), but not the II part (pχ 2 Yates = 0.74); and if it was counted by BJ-algorithm the II part had the biggest amount of n min : by 45.1% more than the both I part and III part (pχ 2 Yates = 0.0002, Table 4). In addition, n max in the I, II, and III part of IFNA17 depending on the each E H cr (from 0.30 to 0.43) were frequently unequal, such as these counted by the decile method (p KWt = 0.0493, p Mt = 0.0764, Figure 3), by the quartile method (p KWt = 0.0335, p Mt = 0.0244, Figure 4), or by BJ-algorithm (p KWt = 0.0220, p Mt = 0.0293, Figure 5); but the base pair amounts of the range "Minimum" dependent on similar conditions in these gene parts had always no differences, e.g., calculated by the decile method (p KWt = 0.491, p Mt = 0.424), by the quartile method (p KWt = 0.386, p Mt = 0.135), or by BJ-algorithm (p KWt = 0.308, p Mt = 0.215).
The difference of the nucleobase quantities of the range "Maximum" between ends and center of the gene was much expressed: pχ 2 Yates was less 0.0001 if n max was counted by the decile method, the quartile method, or BJ-algorithm. Moreover, comparison of the gene ends to the gene center (dependent on the each E H cr ), revealed that the p KWt were less than 0.05 for n max , calculated by the quartile method (0.0475) or BJ-algorithm (0.0295), and p Mt were less than 0.05 for n max , counted by the decile method and BJ-algorithm (both equaled 0.0233). On the other side in this energy range (dependent on the each E H cr ) n min in the ends of IFNA17 comparing to its amount in the center part of the gene had no difference by the decile method (p KWt = 0.505, p Mt = 0.449), by the quartile method (p KWt = 0.382, p Mt = 0.449), or by BJ-algorithm (p KWt = 0.420, p Mt = 0.449).
It should be also noted, that among three gene parts (out of dependent on the energy range), the differences between the total amount of n max and n min were the least in the II part of gene: from 38.3% (range of E H cr : 0.30-0.43, n min > n max ) to 115.1% (range of E H cr : 0.44-0.59, n min > n max ) if counted by the decile method, from 59.8% (range of E H thr : 0.44-0.59, n min > n max ) to 90.5% (range of E H cr : 0.30-0.43, n min > n max )-by the quartile method, or from 7.2% (range of E H cr : 0.30-0.43, n min < n max ) to 11.8% (range of E H cr : 0.44-0.59, n min < n max )-by BJ-algorithm. In contrast of this, the higher differences between the total amount of n max and n min were more often in the III part of gene counted by the decile method, or in the I part-by the quartile method (Table 4).

Discussion
Albeit there should be no difference theoretically in three parts (I, II, and III) of IFNA17 between the amounts of the n min or n max of OSs dependent on the single 2 H/ 1 H replacement, because each part had actually the same number of nitrogenous base pairs in the double-stranded DNA forming it [36], the mathematical modeling demonstrated much more significant differences between their amounts in the described gene parts. The appearance in any gene part of some significant fluctuation of the amount of nucleotide pairs both with the higher P i and the lower P i compared to the P 0 (which correlates to the essential OS frequency) leads to the higher risk of changing the whole DNA functional activity [37]. It is possible due to the fact that a stability of double-stranded DNA supported by the H-bonding between complementary nitrogenous bases and the stacking formed with neighboring bases [38], so the fluctuation opening even of one nucleobase in duplex DNA can involve disruption of the whole stacking in the determined gene part [37,39]. Moreover, decrease of P i in the any gene part lower than essential P 0 level can be also a strong negative factor, because the spontaneous base flipping is the necessary part of some catalyzed processes involving DNA [40,41]. So, for the catalytic activity of DNA glycosylase, 8-oxoguanine DNA glycosylase I, and alkyl adenine DNA glycosylase (which are the important DNA repair enzymes [42]), as well as DNA methyltransferase family enzymes and β-glucosyl transferase (which provide selective modification of DNA [43]) the natural base-flipping frequency needs to be determined. All the mentioned enzymes interact with DNA only according to the nucleobase completely flipped out of the helix structure due to the previous occurrence of OSs in the specific gene region [44][45][46][47]. In addition, the rate of the transcription bubble passages is also dependent on the essential frequency of the open states in the selective part of gene [27,[48][49][50]. Also, both slowdown and rev up of molecular dynamic of DNA due to the change of occurrence OS rate play the significant role in the impairment of DNA-mediated charge transfer (CH-π interaction), which can lead to the dysfunctional activity of repair enzymes [51][52][53], and, consequently, to the increase rate of mutations.
In our study the unequal sensitivity of the different parts of IFNA17 to the single 2 H/ 1 H replacement was demonstrated, which is manifested in the differently expressed changes of their OS frequency. Among three above-described counting approaches of the nucleobase amounts, which can be included in the range "Maximum" or "Minimum" based on the OS frequency occurring after the single 2 H/ 1 H replacement at the specified DNA base pair, the BJ-algorithm should be highlighted. That conclusion can be explained due to the least false positive and false negative results if we are using this algorithm compared to the decile or the quartile methods for the nucleobase counting. For example, in the I and III parts of IFNA17 in the range of E H cr from 0.30 to 0.43, the base pair amounts of "Minimum" calculated according to the decile and quartile methods were different (pχ 2 Yates ≤ 0.0002), but it was because of numerous of the false positive results in data counted by the decile method at E H cr : 0.38 (n min were 13 and 80) and 0.41 (n min were 266 and 308 in the I and III parts respectively, Figure 3), and by the quartile method at E H cr : 0.35 (n min were 283 and 202), 0.36 (n min were 26 and 1), 0.37 (n min were 154 and 36), 0.42 (n min were 66 and 8), and 0.43 (n min were 3 and 25 in the I and III parts respectively, Figure 4). It was observed due to the shifts of P 0 between extreme values of P i : closer to P min or to P max dependent on E H cr (Figures 6 and 7), that produced the false positive and false negative results in the calculations. In contrast to these methods, BJ-algorithm took into account the mentioned shifts of P 0 (Figure 8), that actually avoided the false results and revealed no difference between I and III parts of the gene in this energy range (pχ 2 Yates = 0.954). On the other side, such a new approach cutting down false positive results provides the possibility to find out the significant differences between frequency of OSs in IFNA17 due to the single 2 H/ 1 H replacement in the highlighted gene parts. Although, in the range of E H cr from 0.30 to 0.43 in the I and II gene parts the nucleobase amounts of "Minimum" computed by the decile method had no statistical difference (pχ 2 Yates = 0.85) like the similar nucleobase amounts of the same gene parts calculated by the quartile method (pχ 2 Yates = 0.74), the n min , which were by BJ-algorithm, had reliable difference (pχ 2 Yates = 0.0002), because the last approach diminished total base pair amounts due to the taking away of false positive results: at E H cr 0.38 and 0.41 it was 277 and 209 false n min less than calculated by the decile method ( Figure 3); at E H cr 0.35-0.37 and 0.42-0.43 it was 424 and 347 false n min less than-by the quartile method in I and II parts respectively (Figure 4). However, with that, the new algorithm allowed to decrease false negative results, e.g., when counting n min in the ends of gene and its center part at E H cr 0.32-0.33 and 0.40 (compared to the quartile method, Figures 4 and 5), which increases the reliability of the approach. Another example of the appropriateness of using BJ-algorithm for processing study data is the finding of difference between n min of the ends and center of IFNA17 (pχ 2 Yates < 0.0001), in contrast to the usage of the decile method (pχ 2 Yates = 0.565).   Moreover, more efficient and sustainable usage of BJ-algorithm for identifying n max in the I, II and III parts of IFNA17 confirmed the lowest p-values (p KWt = 0.0220, p Mt = 0.0293) compared to the similar p-values found for n max , which were filtered by the decile and quartile methods (as the pure fractile approach according to the properties of the pure function [54]). Moreover, the p-values of n max between ends and center of the gene were both p KWt and p Mt less than 0.05 (0.0295 and 0.0233, respectively) if they were counted only by BJ-algorithm. The biggest amount of false positive results of n max calculated by the decile method was at E H cr 0.43 (Figure 3), the biggest amount of false negative results of n max filtered by the quartile method was at E H cr 0.36 (Figure 4).
In the range of E H cr from 0.44 to 0.59 the difference among data filtered by three methods was also significant, and BJ-algorithm showed higher efficacy for processing study numbers of OS frequency. For example, the amount of n max in the gene center part filtered by the new algorithm was on 39 n i (which can be considered as false positive results according to the robust estimation) less than n max counted by the quartile method. However, more expressed differences were between n min (Table 4), that explained due to the convergence of P imax and P imin , also P 0 and 0 in the range of higher energies ( Figure 8). For example, the decile method leads to many false positive results of n min at E H cr 0.46 and 0.59, also quartile method-at E H cr 0.55 and 0.59, that provoked among n min of three gene parts (each of which had the same own P i data) completely different p-values: χ 2 B : p > 0.5 (decile method), χ 2 B : p < 0.002 (quartile method), and χ 2 B : p > 0.2 (BJ-algorithm). Contrary to the lower range of E H cr when the fractile methods produced both false positive and false negative results of extreme nucleobase amounts (both n max and n min ), in the higher range of E H cr they much more frequently cause be false positive results in range "Minimum", that stemmed from the following properties of DNA OS frequency: P imax → P imin and P 0 → 0, so all of these p inversely correlated to E H cr ( Table 2). In general, all these practical examples proved the conclusion about the possibility of using BJ-algorithm as the more effective filtering approach of described study data compared to the pure fractile methods. Regarding the above-mentioned information, data will be further considered as filtered only by BJ-algorithm.
Nevertheless, despite the same length of the I, II, and III parts of IFNA17, the study data filtered by each of the methods pointed out to the statistically significant differences of n max and n min of the gene ends and its center, which confirmed heterogeneous sensitivity at the considered gene regions to the single 2 H/ 1 H replacement. It is noteworthy that in the range of E H cr from 0.30 to 0.43 the center part of the gene is the most prone to both fluctuations of n max (more than in 2.73 times) and n min (in 1.45 times) compared to the I and III gene parts. Wherein in the higher energy range (0.44-0.59) II part of IFNA17 significant different from I and III its parts by nucleobase pair amount only in the range "Maximum" (more than in 15.58 times), but there was no difference between n min , because of the convergence extremes of OS frequency in DNA. Note: for each H-bond dissociation energy: the 1 st cross dash is Pimax, the 2 d cross dash is bottom of the range "Maximum", the 3 d cross dash is top of the range "Minimum"; the 4 th cross dash is Pimin; red line is OS occurrence frequency under condition when all hydrogen bonds in DNA are 1 H (P0).
All described fluctuations of nmax and nmin after single 2 H/ 1 H replacement can result in DNA damage due to the slowdown in the rate of bubble forming (especially non-superhelical stressinduced denaturation bubbles having more than 4 base-pairs), including the promoter regions and binding points of specific proteins (e.g., DNA repair enzymes [41,55]), or because of the overwhelming rate of flipping-out DNA nucleobases that can be accompanied by increasing both the speed of modification of their chemical structure and mismatched protein-DNA interaction [56][57][58][59].
Our data clearly demonstrate the possibility of higher risk of gene lesion in all its length due to single 2 H/ 1 H replacement predominantly in the center part of IFNA17 in all studied critical energies (from 0.30 to 0.59), which show abnormal rate of occurring OSs, that manifested more often its increasing (by 24%) compared to decreasing. According to the occurrence rate of Pimax in the range of E H cr from 0.44 to 0.59, the single deuterium introduction in the III part of gene actually did not cause any changes compared to the natural frequency of OSs as opposed to the single 2 H/ 1 H replacement in I or II parts. So, after deuterium introduction in the III part, predominated nmin, which can lead to the sharp slowdown of DNA bubble forming, which has more than 4 nucleobases participating in generation of the preexisting state of denaturation, that is required by specific DNA-binding enzymes [56,59], and can be a predictor of function lesion of DNA repair system. It was well known that different DNA loci have inequality in the adenine-thymine/guaninecytosine (A-T/G-C) ratio. In IFNA17 I part has 163 A-T and 164 G-C base pairs, II part has 186 A-T and 140 G-C base pairs, and III part has 236 A-T and 91 G-C base pairs [30]. Despite the fact, that I part has the minimum number of A-T pairs and the III part has the maximum number of A-T pairs (compared to each other), the most significant change of OS frequency in IFNA17 was found out after single 2 H/ 1 H replacement in the II part of gene (which contained intermediate AT percentage) in the all studied energy ranges (it was observed under both E H cr diapasons: from 0.30 to 0.43 and from 0.44 . Note: for each H-bond dissociation energy: the 1st cross dash is P imax , the 2nd cross dash is bottom of the range "Maximum", the 3rd cross dash is top of the range "Minimum"; the 4th cross dash is P imin ; red line is OS occurrence frequency under condition when all hydrogen bonds in DNA are 1 H (P 0 ).
All described fluctuations of n max and n min after single 2 H/ 1 H replacement can result in DNA damage due to the slowdown in the rate of bubble forming (especially non-superhelical stress-induced denaturation bubbles having more than 4 base-pairs), including the promoter regions and binding points of specific proteins (e.g., DNA repair enzymes [41,55]), or because of the overwhelming rate of flipping-out DNA nucleobases that can be accompanied by increasing both the speed of modification of their chemical structure and mismatched protein-DNA interaction [56][57][58][59].
Our data clearly demonstrate the possibility of higher risk of gene lesion in all its length due to single 2 H/ 1 H replacement predominantly in the center part of IFNA17 in all studied critical energies (from 0.30 to 0.59), which show abnormal rate of occurring OSs, that manifested more often its increasing (by 24%) compared to decreasing. According to the occurrence rate of P imax in the range of E H cr from 0.44 to 0.59, the single deuterium introduction in the III part of gene actually did not cause any changes compared to the natural frequency of OSs as opposed to the single 2 H/ 1 H replacement in I or II parts. So, after deuterium introduction in the III part, predominated n min , which can lead to the sharp slowdown of DNA bubble forming, which has more than 4 nucleobases participating in generation of the preexisting state of denaturation, that is required by specific DNA-binding enzymes [56,59], and can be a predictor of function lesion of DNA repair system. It was well known that different DNA loci have inequality in the adenine-thymine/guanine-cytosine (A-T/G-C) ratio. In IFNA17 I part has 163 A-T and 164 G-C base pairs, II part has 186 A-T and 140 G-C base pairs, and III part has 236 A-T and 91 G-C base pairs [30]. Despite the fact, that I part has the minimum number of A-T pairs and the III part has the maximum number of A-T pairs (compared to each other), the most significant change of OS frequency in IFNA17 was found out after single 2 H/ 1 H replacement in the II part of gene (which contained intermediate AT percentage) in the all studied energy ranges (it was observed under both E H cr diapasons: from 0.30 to 0.43 and from 0.44 to 0.59, Table 4). All of this definitely points out that not only total number (or percentage) of AT nucleobases influences the OS frequency [60], but the specific sequence of A-T and G-C bases, or DNA pattern, which can lead to the both increase and decrease of OS amount in the gene. Thus, all of described cases of dynamics of the double-stranded DNA confirmed the heterogeneous sensitivity of IFNA17 to the 2 H/ 1 H exchange dependent on an energy of external influences on the gene. The study also showed that the single deuterium introduction in IFNA17 even in the III gene part, which is out of the coding region (from 50th to 619th nucleobases [61,62]), can influence the transcription speed by modifying the dynamics of a double-stranded DNA due to the decreasing of OS amounts in other gene parts. Moreover, the coding region of IFNA17, which is more relevant to I and II parts, is the most sensitive part of the gene by the single 2 H/ 1 H replacement. However, our study had some limitations. One of them was that all results have been obtained in the frameworks of the mathematical model based on the Newton's equations and representative of the Cauchy problem for the system of 2n ordinary differential equations, which do not entirely take into account the contribution of DNA-protein interactions. The second limitation is increasing data variability in the low energy range (less than 0.30) due to the properties of the system of 2n ordinary differential equations. Nevertheless, we do not consider these limitations as an insurmountable obstacle, and think that the approach can be generalized and used to more complex models of DNA dynamics.

Conclusions
In general, the mathematical modeling study demonstrated significant inequality (dependent on single 2 H/ 1 H replacement in DNA) among three parts of gene similar in length of the frequency of occurrence of the OSs both with the higher P i and the lower P i compared to the natural P 0 . The bases pairs with P i significantly higher than P 0 can provide the bigger risk of mutation (due to a change in the ratio of the equilibrium constants of the denaturation bubble opening reaction in various DNA parts, including the promoter regions, and the thermodynamic indices at binding points of specific enzymes [27,63,64]), whereas the nucleotide pairs with Pi significantly lower than P 0 can lead to the harsh slowdown of transcription or replication compared to the natural level of failures in those essential processes [49,50,55]. In addition, in the study it was found out that through the range of H-bond dissociation energy more than 0.58 and approximated to 0.59 the DNA OSs were not occurred in IFNA17, and only close states of nitrogenous base pairs were revealed in the whole gene.
In this paper the new convenient approach of the analysis of the abnormal frequency of OSs in the different part of IFNA17 was presented, which took into account both rising and decreasing of them that allows to make a prediction of the functional instability of the specific DNA regions. One advantage of the new algorithm is diminishing the number of both false positive and false negative results in data filtered by this approach compared to the pure fractile methods, such as deciles or quartiles. Thus, study data pointed out to the statistically significant differences of nucleobase pair amounts, which can be included in the extreme ranges ("Maximum" or "Minimum") based on the frequency of OS occurrence, in the gene ends and its center, so the dynamics of a double-stranded DNA confirmed its heterogeneous sensitivity to the single 2 H/ 1 H replacement dependent on energy of external influences on the gene. It was also shown that the single deuterium exchange in IFNA17 even out of the coding region (from 50th to 619th base pairs), e.g., in the III gene part, can probably influence the transcription speed by modifying the dynamics of a double-stranded DNA due to the decreasing of OS amounts in other gene parts. . . . . . .