Influence of Single Deuterium Replacement on Frequency of Hydrogen Bond Dissociation in IFNA17 under the Highest Critical Energy Range

The effect of single substitutions of protium for deuterium in hydrogen bonds between pairs of nitrogenous bases on the open states occurrence probability at high critical breaking energies of these bonds has been studied. The study was carried out using numerical methods based on the angular mathematical model of DNA. The IFNA17 gene was divided into three approximately equal parts. A comparison of the open states occurrence probability in these parts of the gene was done. To improve the accuracy of the results, a special data processing algorithm was developed. The developed methods have shown their suitability for taking into account the occurrence of open states in the entire range of high critical energies. It has been established that single 2H/1H substitutions in certain nitrogenous bases can be a mechanism for maintaining the vital activity of IFNA17 under critical conditions. In general, the developed method of the mathematical modeling provide unprecedented insight into the DNA behavior under the highest critical energy range, which greatly expands scientific understanding of nucleobases interaction.


Introduction
Experimental study of the DNA structure is difficult because of several problems; the most important of which is the limitation of the spatial resolution of available research tools [1,2]. Despite the development of methods for studying single molecules, they have a few limitations for studying the mechanics of DNA [3,4]. For this reason, mathematical modeling is one of the main modern methods for studying DNA molecular dynamics, its mechanical movements, open states, denaturation bubbles [5][6][7]. This method, despite several simplifications, allows us to consider various aspects of DNA functioning with great accuracy [8][9][10]. In the article [11], using a combination of atomic force microscopy (AFM) and modeling of atomic molecular dynamics (MD), an attempt was made to experimentally observe the dynamics of open states. The experiment was carried out on short ring DNA, while many microphotographs were taken, after which, using a custom Python script, the images were corrected for further morphological analysis. As a result, denaturation bubbles were recorded, which, in turn, leads to DNA compaction. This indicates the importance of combining the experimental and mathematical modeling methods, which allows us to obtain much more information about the dynamics and open states of the DNA molecule [12].
Earlier, we found that the ingress of a deuterium atom into hydrogen bonds between pairs of nitrogenous bases increases the probability of occurrence of open states by 0.22-0.60% [13]. In addition, it has been shown that the probability of the bubbles formation (length from 12 to 27 nucleotides) depends on the localization of the deuterium atom in the DNA molecule and may differ significantly from the probability of the occurrence of open states in general [14]. The participation of deuterium atoms in the formation of hydrogen bonds in the DNA molecule can cause a change in the time of transmission of genetic information, thus, even a slight change in the isotopic state of the medium can affect changes in metabolic processes in living systems [15][16][17]. It is known that non-radioactive isotopes of biogenic elements ( 2 H/ 1 H, 13 C/ 12 C, 15 N/ 14 N, 18 O/ 17 O/ 16 O) have a significant effect on the rate of biochemical reactions, physiological processes, growth, and development of unicellular and multicellular living organisms with different levels of organization of energy metabolism and metabolic rate [18][19][20][21][22][23][24][25].
The aim of this work is to study the effect of a single deuterium substitution at the frequency of hydrogen bond dissociation in IFNA17 in the range of the highest critical energies, based on the mechanical model of DNA [26], without simplifications and averaging [27,28].

Results
On Figure 1 the graphs of the angular deviations of the 1-st chain of DNA molecule nitrogenous bases under the first E H cr , which equals 0.581 × 10 −22 N·m over period from 0 to 3.0 × 10 −10 s are presented.
(AFM) and modeling of atomic molecular dynamics (MD), an attempt was made to experimentally observe the dynamics of open states. The experiment was carried out on short ring DNA, while many microphotographs were taken, after which, using a custom Python script, the images were corrected for further morphological analysis. As a result, denaturation bubbles were recorded, which, in turn, leads to DNA compaction. This indicates the importance of combining the experimental and mathematical modeling methods, which allows us to obtain much more information about the dynamics and open states of the DNA molecule [12].
Earlier, we found that the ingress of a deuterium atom into hydrogen bonds between pairs of nitrogenous bases increases the probability of occurrence of open states by 0.22-0.60% [13]. In addition, it has been shown that the probability of the bubbles formation (length from 12 to 27 nucleotides) depends on the localization of the deuterium atom in the DNA molecule and may differ significantly from the probability of the occurrence of open states in general [14]. The participation of deuterium atoms in the formation of hydrogen bonds in the DNA molecule can cause a change in the time of transmission of genetic information, thus, even a slight change in the isotopic state of the medium can affect changes in metabolic processes in living systems [15][16][17]. It is known that non-radioactive isotopes of biogenic elements ( 2 H/ 1 H, 13 C/ 12 C, 15 N/ 14 N, 18 O/ 17 O/ 16 O) have a significant effect on the rate of biochemical reactions, physiological processes, growth, and development of unicellular and multicellular living organisms with different levels of organization of energy metabolism and metabolic rate [18][19][20][21][22][23][24][25].
The aim of this work is to study the effect of a single deuterium substitution at the frequency of hydrogen bond dissociation in IFNA17 in the range of the highest critical energies, based on the mechanical model of DNA [26], without simplifications and averaging [27,28].

Results
On Figure 1 the graphs of the angular deviations of the 1-st chain of DNA molecule nitrogenous bases under the first E , which equals 0.581 × 10 −22 N·m over period from 0 to 3.0 × 10 −10 s are presented.  The open state (OS) occurrence frequency was counted by using above-mentioned and modified Basov-Jimack algorithm, according to which, for the highest critical energy range (from 0.581 × 10 −22 N·m to 0.589 × 10 −22 N·m), on Figure 2 the dynamics of the OS occurrence in studied gene dependent on the H-bond dissociation energy under natural condition and after the single 2 H/ 1 H replacement are showed. This data, more than anything else, demonstrates that, in the highest critical energy range, the OS occurrence frequency under natural condition, when all hydrogen bonds in DNA nucleotides are 1 H, is always lower than frequency of OS occurrence in the range "Maximum", when is the single Figure 2. Dynamics of the open state (OS) occurrence in gene encoding interferon alpha 17 dependent on the H-bond dissociation energy under natural condition and after the single 2 H/ 1 H replacement (with gradation of OS occurrence frequency by modified Basov-Jimack algorithm). Note: for each H-bond dissociation energy: the 1st cross dash is , the 2nd cross dash is the bottom of the range "Maximum", the 3d cross dash is the top of the range "Minimum", which was calculated by using non-modified Basov-Jimack algorithm [29]; the red line is OS occurrence frequency under condition when all hydrogen bonds in DNA are 1 H (P 0 ). . Note: for each H-bond dissociation energy: the 1st cross dash is P i max , the 2nd cross dash is the bottom of the range "Maximum", the 3d cross dash is the top of the range "Minimum", which was calculated by using non-modified Basov-Jimack algorithm [29]; the red line is OS occurrence frequency under condition when all hydrogen bonds in DNA are 1 H (P 0 ).

Figure 3.
Distribution of nucleobase pairs (which was counted by modified Basov-Jimack algorithm in the different parts of IFNA17), leading after the single 2 H/ 1 H replacement to the extreme frequencies of OS and CSNB occurrences. Note: red dot is the location of the deuterium atom in the DNA molecule, which leads to the maximum probability of OS occurrence (range); green dot is the location of the deuterium atom in the DNA molecule, which leads to the CSNB occurrence (range).

Figure 3.
Distribution of nucleobase pairs (which was counted by modified Basov-Jimack algorithm in the different parts of IFNA17), leading after the single 2 H/ 1 H replacement to the extreme frequencies of OS and CSNB occurrences. Note: red dot is the location of the deuterium atom in the DNA molecule, which leads to the maximum probability of OS occurrence (range); green dot is the location of the deuterium atom in the DNA molecule, which leads to the CSNB occurrence (range).  I  163  164  II  186  140  III  236  91 Note: CSNB is closed states of nitrogenous bases (PCS = 0), A-adenine, T-thymine, G-guanine, C-cytosine.
On Figure 3, in the different parts of IFNA17, showed the distribution of nucleobase pairs leading to the maximum probability of OS and CSNB occurrences after the single 2 H/ 1 H replacement under energy diapason (E H cr ) from 0.581 × 10 −22 N·m to 0.589 × 10 −22 N·m (henceforth 0.581-0.589).
In our previous study for 0.30-0.58 × 10 −22 N·m energy range [29], the inequality in the frequency of the open states occurrence because of the single 2 H/ 1 H replacement in gene encoding interferon alpha 17 (IFNA17, which consists of 980 nucleotide pairs [30]) was proved through it conditionally divided into three equal parts: from the 1st to the 327th nucleobases (I part), from the 328th to the 653th nucleobases (II part) and from the 654th to the 980th nucleobases (III part). This approach allowed us to contemplate the influence on probabilities (P i ) of open state occurrences between different nitrogenous bases in double-stranded DNA dependent on the single 2 H/ 1 H replacement in basepair of each gene region. Although these studied parts of IFNA17 had approximately the same number of base pairs, the A-T/G-C ratios were significantly different in each part [31], which are presented in more detail in Table 2. Table 2. Distribution of adenine-thymine (A-T) and guanine-cytosine (G-C) base pairs in different parts of IFNA17. According to this fact, for the range of the energies from 0.581 × 10 −22 N·m to 0.589 × 10 −22 N·m both quantity of closed states of nitrogenous bases (CSNB) and number of nucletides of the range "Maximum" (n max ) for different nitrogenous bases in each of three parts of IFNA17 under single 2 H/ 1 H replacement via BJ-algorithm were calculated, and obtained data are presented in Table 1 for CSNB and Table 3 for n max . Table 3. Quantity of A-T and G-C nucleobases, included in the range "Maximum", in each of three parts of IFNA17 dependent on single 2 H/ 1 H replacement.  It is worth noting that the change in energy, which was taken into account when calculating CSNB (Table 1) and n max (Table 3) in the highest energy range, always equaled 0.001.

Part of IFNA17 Nucleobase Pair Quantity A-T Ratio (%) G-C Ratio (%)
In our study was found out that quantity of CSNB and n max of OSs for the each gene part were different in the above-described energies (E H cr ) and depend on A-T/G-C ratio (Tables 1 and 3). In the range of E H cr from 0.581 to 0.589 the highest total number of n max counted by the BJ-algorithm was when the single 2 H/ 1 H replacement had taken place in the II part of IFNA17. The number of i, which were included in the range "Maximum" in this gene part was in 3.21 times higher than in the I part (χ 2 B : p < 0.001); whereas counted n max was equal to 0.0 for any of studied energies when 2 H/ 1 H replacing had occurred in the III part ( Figure 3). Albeit the A-T nitrogenous bases prevalence over G-C ones in IFNA17 (in 1.48 times), the distribution of the total n max in the whole gene due to the isotopic substitution at the A-T and G-C nucleotides showed the predominance of 2 H/ 1 H replacement at the cytosine and guanine bases, the number of which was in 3.67 times higher, than the number of 2 H-substituted A-T (pχ 2 B < 0.0001). Wherein, the difference between n max of 2 H/1H-substituted A-T and G-C nucleotides was much more pronounced in the gene I part (in 8.5 times, Table 3) compared to II part (in 3.03 times, Table 3). In addition, the overwhelming majority of 2 H/ 1 H-replacement at the G-C pairs, leading to the n max , arose in the diapason of E H cr from 0.581 to 0.585, so that the number of them was higher in 7.56 times than the nmax via the single 2 H-substituted adenine or thymine in the same range of energies (pχ 2 B < 0.0001). For E H cr from 0.586 to 0.589, there is a less pronounced difference between deuterium substituted A-T and G-C n max , although the number of the last is more in 2.29 times (pχ 2 B < 0.00001, Table 3). The calculated quantity of CSNB firstly turned higher in the I part of IFNA17 (under E H cr from 0.581 to 0.582) compared to the II and III parts, which having nothing of these ( Figure 3). Further, under E H cr from 0.583 to 0.584 the CSNBs arise in the gene II part; but even at energy equal to 0.584 their number keeps less in 2.75 times than in the I part. Only upon reaching the influence of 0.585 energy in the III part appears some CSNBs, the number of which was less compared to the I and II parts in 1.42 and 1.69 times, respectively (χ 2 B : p < 0.03, Figure 3). However, starting with E H cr equal to 0.586, the quantity of CSNBs in the III part was grown very abruptly, reaching 100% of the nucleotides of this gene part, and it exceeded the CSNB numbers, which were almost equal to each other, in the I and II parts of IFNA17 in 1.51 and 1.56 times, respectively (χ 2 B : p < 0.006, Figure 3). In addition, as in the coming into being of n max , the influence of A-T and the G-C ratio was explicit for generation of CSNBs under the single deuterium substitution in each part of the studied gene. In the 0.581-0.582 range of critical energies, the CSNBs were initiated and grew up specific via the 2 H-substituted G-C pairs, and, moreover, in the 0.583-0.584 E H cr range the 2 H-substituted G-C outnumbered the 2 H-substituted A-T nitrogenous bases and were accounted at least 67% of the total nucleobases leading to the CSNB occurrence in the whole gene (pχ 2 B = 0.0056, Table 1). However, at the E H cr equal to 0.585 this difference of CNBS numbers between 2 H-substituted A-T and G-C was leveled (pχ 2 B = 0.977). It should be noted, that, under 0.589 E H cr , CSNB making was completed for the each nucleobase pairs of the whole IFNA17, and it observed when six G-C pairs had turned into close states in the gene II part specifically (pχ 2 B = 0.0108). In addition, it was found that than the lower A-T/G-C ratio in the gene part, the CSNBs arose there earlier and under lower critical energies. So, for example, in the range of E H cr from 0.581 to 0.585 the Spearman correlation coefficient between A-T/G-C ratio and CSNB numbers was −0.547 (p = 0.035, Figure 3). The fractions of the 2 H-substituted G-C bases in the initiative numbers of CSNB in the each gene parts were following: 100% in the I part under E H cr equal to 0.581 (pχ 2 B = 0.2547), 100% in the II part under E H cr equal to 0.583 (pχ 2 B = 0.0359), and only 30% in the III part under E H cr equal to 0.585 (pχ 2 B = 0.771, Table 1), that could be because of the last gene part had the least number of the G-C nucleobases ( Table 2). For the highest acceleration rates of the CSNB numbers, with reaching them approximately half of the total nucleobases in each gene part (but less than all of the nitrogenous bases in the each part) were provided more often due to the 2 Table 1).

Discussion
The influence of the isotopic 2 H-effect on nucleic acids structure and their functions is being actively studied, and it is relevant for various fields of science [19,[32][33][34][35][36][37][38][39][40]. So, our data allows to understand the mechanism of the possible changes in the OS occurrence frequency in IFNA17 by the single 2 H/ 1 H-replacement under specific energy range more precisely. It is well known about deuterium (D) ability to alter structure and activity of DNA and RNA via changing base-pairing interaction strength by the difference in the energy, stability and geometry between H-and D-bonds, and also through the inequality between 1 H and 2 H in the hydrogen bond zero-point vibration energies [32,35,36,[41][42][43][44][45][46]. In the present study one of the most vital processes for DNA functioning, which can depend on isotopic 2 H/ 1 H-exchange, was considered-generation of OSs, surplus or deficiency of which can provoke various impairments of DNA stability. It is worth noting that occurrence of nucleotide open states in DNA is very critical conformational transition, which is necessarily needed for both hydrogen exchange and intermolecular interaction of proteins with specific DNA target, which, consequently, provide the realization of the vital function of nucleic acids [47]. In addition, expressing concisely, CSNB are logically inappropriate for the specific implementation of some essential processes in cell due to the formation in DNA many excessively unbreakable sections, which totally prevent the hydrogen exchange between proteins and DNA. Additionally, those energies, under which gene has one or more CSNB (which can equal in summary up to the entirety number of its nucleotides), form the highest critical energy range for the certain gene. Therefore, the study of the influence of 2 H/ 1 H exchange on the molecular dynamic of DNA for energy diapason (or the highest energy range), which initials the rising of the nucleotide closed states in IFNA17, and also taking into consideration A-T/G-C ratio in its different parts, is the great interest for a deeper understanding of the biological functioning of nucleic acids under critical environmental effects [41]. So, according to the data obtained for the studied IFNA17 [29], the highest energy range equals to E H cr from 0.581 to 0.589 for it (Table 1, Figure 3), and that was contemplated in this work.
According to the developed mathematical modeling it was found in the study, that the more increased stability of IFNA17 observed in its parts rich in G-C base pairs, which can be associated with the specific DNA folding of these regions providing more stabilizing energetic contributions and earlier transition to closed states of nucleobases under lower critical energies (starting from the E H cr equals to 0.581 in the gene I part). This statement is confirmed by the significant and negative Spearman correlation coefficient between A-T/G-C ratio and CSNB numbers in the each of the three studied parts of IFNA17. It is well known that even the change of closing nitrogenous bases from G-C to C-G can stabilize DNA hairpin structure approximately 2 kcal/mol and lead to the significant enhance of its melting temperature [48,49]. The similar phenomena are most possible associated with the formation in the alternative DNA structures additional H-bonds, for instance, due to the slightly twisted some nitrogenous base pairs. In addition, it is possible through the occurrence of electrostatically favorable contacts between the partially positively charged amino groups and partially negative oxygen-containing groups with the subsequent formation of different types of loops or other energetically auspicious DNA structures [50,51]. Another but less influenceable reason of the arising CSNB in the gene is increasing of van der Waals packing forces, which changed at the specific conformations of DNA. In opposing to CSNB, the increase OS occurrence frequency, which stems from the presence of the single 2 H/ 1 H-substituted nucleotide in the specific gene region, requires weakening of hydrogen bonds in DNA, reduction of electrostatic attraction between charged chemical groups and decreasing of van der Waals packing forces, that occurred firstly and more often when 2 H/ 1 H replacement takes place at the cytosine and guanine bases, especially in the II part of IFNA17. It can be induced by different reasons: some conformational change of the gene, specific molecular folding, and reducing rate of the rare tautomer formation within DNA because of the isotopic 2 H/ 1 H exchange in nitrogenous bases [52,53]. In our study was showed, that the single 2 H/ 1 H replacement at G-C pairs has the most influence on the initial arise both n max and n CSNB . Moreover, the last OSs of IFNA17 under E H cr = 0.588 × 10 −22 N·m were presented only G-C pairs ( Table 1). The peculiarity of the gene III part with a significant higher content of A-T is the concise turning OSs into CSBNs within two critical energies, such as a collapse (Figure 3). It can point out the higher destabilizing of DNA region enriched A-T nucleobases and risk of its sharp dysfunction under critical conditions [54], nevertheless, the initial OS frequency changes occur foremost because of the 2 H-substitution at G-C pairs.
In addition, the above modification of the BJ-algorithm was due to the fact that, in the highest critical energy range, there was a progressive accumulation of the number of CSNB (that was starting from E H cr eaqualed 0.581 when for the first time the minimum P i limit was reached: P i = 0.0). Further, this did not allow the differential consideration of the contribution of the de novo emerging CSNBs (under non-modified BJ-algorithm [29]), which were beginning continually rising from E H cr equaled 0.582 ( Figure 3). So that, since the subsequent CSNBs additively accumulated in the total number of the previously formed CSNBs, the role of each emerging CSNB has been indistinguishable in the total number of CSNBs ( Figure 3).
Bearing in mind following prerequisites: the contribution of earlier arisen CSNB is clearly higher in failure of the molecular dynamics of the gene, than the later ones, and for the newly formed CSNB it is also inversely proportional to the summary number of the ones in gene, as well as the number of CSNBs are limited by the total number of gene nucleobase pairs the modified BJ-algorithm can required additional below-presented changes for this study: i range "Minimum" (CSNB): if 0 < n g CSNB /n g ≤ 1: where: CSNB is closed state of nitrogenous bases: PCSNB = 0; n x CSNB is number of closed states of nitrogenous bases in the gene part x under certain E H cr ; n g CSNB is total number of closed states of nitrogenous bases in the whole gene under certain E H cr ; n g is number of nitrogenous base pairs in the whole gene; If P i ≤ P i min + 3 4 P 0 − P i min : i range from Q 2 -min to Q 4 -min (Q 2 -Q 4 -min) [29]. The new approach can not only take into account the bigger role earlier formed CSNBs for the certain gene functioning (due to accounting for n x CSNB E H cr depending on n g CSNB E H cr /n g , but also avoid the cumulative effect distorting the total CSNB number (especially for the higher energies: E H cr ≥ 0.586, Figures 3 and 4) by using square function for the component: 1 − n g CSNB E H cr /n g . Foremost this approach allows us to reduce the number of false positive data when evaluating CSNBs (Figures 3 and 5). Other benefits of the changes in the new approach of BJ-algorithm and some disadvantages of decile and quartile methods for selection of the nucleobases for the maximum and minimum ranges were presented in details in our earlier work [29].
For instance, if we confront statistics of the data, processed by non-modified BJalgorithm and the new approach, it can be found when comparing their results, there is no significant difference between the sums of n CSNB , which were counted by the nonmodified BJ-algorithm [29], in the second (II) and third (III) parts of the gene throughout 0.581-0.589 energy range (n (part-II) = 1371 bases, n (part-III) = 1415 bases, pχ 2 Yates = 0.3927), but there is significant difference on 205.9% between the sums of n CSNB , higher in the II part compared to III part, are calculated by the new approach in the same energy diapason (n (part-II) = 105 bases, n (part-III) = 51 bases, pχ 2 Yates < 0.0001). It was reached due to the strong reduction of the false positive results because of the above-mentioned new approach of the modified BJ-algorithm (Figures 3 and 5). Additionally, when CSNB numbers were counted by the new approach of the modified BJ-algorithm in the whole gene and for each E H cr throughout the energy diapason from 0.581 to 0.588, the strong correlation between the total n CSNB was revealed (because of 2 H-substitution at both A-T and G-C) and n CSNB due to the 2 H-substitution only at G-C nucleobases, with Spearman's coefficient equaled to 0.881 (p = 0.00385), which proved the far higher ability of 2 H-substituted G-C produced more evenly and incessantly the closed states in IFNA17 compared to 2 H-substituted A-T bases ( Figure 6). the component: 1 − CSNB (E )/ . Foremost this approach allows us to reduce the number of false positive data when evaluating CSNBs (Figures 3 and 5). Other benefits of the changes in the new approach of BJ-algorithm and some disadvantages of decile and quartile methods for selection of the nucleobases for the maximum and minimum ranges were presented in details in our earlier work [29]. Figure 4. Fraction of nucleobase pairs, which was counted by modified Basov-Jimack algorithm in the different parts of IFNA17, leading after the single 2 H/ 1 H replacement to the CSNB occurrences under the highest critical energy range. Note: fraction of CSNB was calculated as CSNB number under specific E divided by the total nucleobase number in the each part of gene; fraction of CSNB in the gene was calculated as CSNB number in the whole gene under specific E divided by the total nucleobase number in the gene (n = 980); part I is from 1st nucleotide to 327th nucleotide; part II is from 328th nucleotide to 653d nucleotide; part III is from 654th nucleotide to 980th nucleotide.
For instance, if we confront statistics of the data, processed by non-modified BJ-algorithm and the new approach, it can be found when comparing their results, there is no significant difference between the sums of CSNB , which were counted by the non-modified BJ-algorithm [29], in the second (II) and third (III) parts of the gene throughout 0.581-0.589 energy range (n(part-II) = 1371 bases, n(part-III) = 1415 bases, 2 = 0.3927), but there is significant difference on 205.9% between the sums of CSNB , higher in the II part compared to III part, are calculated by the new approach in the same energy diapason (n(part-II) = 105 bases, n(part-III) = 51 bases, 2 < 0.0001). It was reached due to the strong reduction of the false positive results because of the above-mentioned new approach of the modified BJ-algorithm (Figures 3 and 5). Additionally, when CSNB numbers were counted by the new approach of the modified BJ-algorithm in the whole gene and for each E throughout the energy diapason from 0.581 to 0.588, the strong correlation between the total CSNB was revealed (because of 2 H-substitution at both A-T and G-C) and CSNB due to the 2 H-substitution only at G-C nucleobases, with Spearman's coefficient equaled to 0.881 ( = 0.00385), which proved the far higher ability of 2 H-substituted G-C produced more evenly and incessantly the closed states in IFNA17 compared to 2 H-substituted A-T bases ( Figure 6).   Figure 5. Number of nucleobase pairs, which was counted by the new approach of modified Basov-Jimack algorithm in the different parts of IFNA17, leading after the single 2 H/ 1 H replacement to the CSNB occurrences under the highest critical energy range. Note: fraction of CSNB in the gene was calculated as CSNB number in the whole gene under specific E divided by the total nucleobase number in the gene (n = 980) and then multiplying this ratio by 100%. In general, the developed method of the mathematical modeling provide unprecedented insight into the DNA behavior under the highest critical energy range, which greatly expands scientific understanding of nucleobases interaction [50]. Under the highest critical energy range, these described fluctuations of CSNB after single 2 H/ 1 H-substitutions in the different gene parts can provoke DNA dysfunction, for example, because of the sharp slowdown in the rate of the generation of vital non-superhelical denaturation bubbles, which have more than four nucleobase-pairs in the promoter regions and specific In general, the developed method of the mathematical modeling provide unprecedented insight into the DNA behavior under the highest critical energy range, which greatly expands scientific understanding of nucleobases interaction [50]. Under the highest critical energy range, these described fluctuations of n CSNB after single 2 H/ 1 H-substitutions in the different gene parts can provoke DNA dysfunction, for example, because of the sharp slowdown in the rate of the generation of vital non-superhelical denaturation bubbles, which have more than four nucleobase-pairs in the promoter regions and specific binding points for DNA repair proteins [55,56] that can be a predictor of DNA repair system lesion. In addition, obtained data demonstrate the possibility of higher risk, which increase abruptly in the whole energy range from 0.581 to 0.587 of altering and impairment of IFNA17 due to the single 2 H/ 1 H exchange in the A-T richest part compared to the both I or II parts, especially under E H cr more than 0.585 (Figures 3 and 4). It should be noted, that in opposition to the natural frequency of OSs in the gene under the energy higher 0.586, after the single 2 H-substitution at G-C in the II part observed lower n CSNB , so this can provide more probability of DNA bubble occurrence, be one of nucleic acid adaptation mechanisms to critical conditions, and matter to the evolution process of living system [57]. According to the developed mathematical modeling and modified BJ-algorithm it was found, that, under the E H cr from 0.586 to 0.588 the occurrence P i rates, which matches the range "Maximum", due to the 2 H/ 1 H exchange at G-C is always significant more than ones due to the same isotopic exchange at A-T, that can consider as the positive participating of the 2 H-substituted guanine and cytosine in nucleic acid adaptation in opposing to the 2 H-substituted adenine and thymine, which give no advantages to DNA denaturation process. For example, it can influence the rate of specific DNA bubble forming, which participating in generation of the pre-existing state of denaturation, that is required by specific DNA-binding enzymes [58,59]. The relationship study between total n max and n max produced due to the 2 H/ 1 H substitution at G-C in the I and II parts of IFNA17, which had more than one 2 H-substituted bases from range "Maximum" in opposing to III part, for each E H cr throughout the energy diapason from 0.581 to 0.588 revealed Spearman's coefficient equals to 0.994 (p < 0.00001), that pointed out the obvious dominant role of 2 H-substituted G-C bases in the generation n max in these gene parts compared to A-T ones ( Figure 7). Moreover, the study demonstrated, that the single deuterium introduction in the gene III part, which is out of the coding region from 50th to 619th base pairs [60,61], had the most close data of the OS occurrence rate compared to the data of OS occurrence frequency under condition when all hydrogen bonds in IFNA17 are 1 H, that was showed clearly the least possibility of 2 H/ 1 H replacement in this part to influence on gene adaptation process under critical conditions positively. So, below are the main inferences of our work, confirming in detail the above-mentioned conclusion: 1. In some cases, the single 2 H/ 1 H replacement result in positive value of OS occurrence frequency in the IFNA17 throughout E diapason from 0.586 × 10 −22 N·m to 0.588 × 10 −22 N·m opposing to the OS occurrence frequency, when all hydrogen bonds in DNA nucleotides are 1 H, which is always equaled to 0.0 after E exceed 0.585 × 10 −22 N·m. This underlines that at least 22.7% of the total number 2 H-substituted nucleobases can reduce molecular interaction in the studied gene and increase the hydrogen bond dissociation, foremost in its part from 328 to 653 nitrogenous bases; 2. The counted occurrences of the OSs were much higher when the single 2 H/ 1 H replacement had taken place at nucleobases of the middle part IFNA17 (from 328 to 653 nucleotides) compared to its other parts, and this had a strong prevalence rate in the G-C pairs; 3. The lowest rate of the OS occurrence was under the single deuterium substitution at the nitrogenous bases in the gene III part (from 654 to 980 nucleotides), which was also too rich in A-T pairs (72.2%) compared to the other parts of IFNA17, so that the calculated was equal to 0.0 for all of the studied critical energies (from 0.581 × 10 −22 N·m to 0.589 × 10 −22 N·m); 4. Sum of was less significant when the single 2 H/ 1 H replacement occurred at the A-T nucleobase pairs compared to the G-C ones in the I and II parts of IFNA17, and the relationship between total and at G-C in these parts for each E throughout the energy diapason from 0.581 × 10 −22 N·m to 0.588 × 10 −22 N·m was strong ( Spearman = 0.994) that proves the obvious dominant role of 2 H-substitutied G-C bases in the generation compared to the A-T; 5. Earliest CSNBs (n = 3) arose under E equal to 0.581 × 10 −22 N·m when IFNA17 2 Figure 7. Numbers of nucleobases from the range "Maximum", which were counted by the modified Basov-Jimack algorithm, in the I and II parts of IFNA17 throughout 0.581-0.588 energy diapason. Note: Total n max in each gene part was calculated as n max number due to the single 2 H/ 1 H-exchange both at A-T and G-C pairs in total; G-C n max in each gene part was calculated as n max number due to the single 2 H/ 1 H-exchange only at G-C pairs.
Above-described cases of the double-stranded DNA dynamics obviously confirmed that under the highest critical energy diapason, the heterogeneous sensitivity of IFNA17 to isotopic 2 H/ 1 H exchange took place and frequency of OS occurrence strongly depended on not only the single 2 H-modified nucleobase in the certain gene region (I, II, or III) but, additionally, on the A-T/G-C-ratio in each of its part; the last was, especially, correlated to the quantity of initial CSNB, that was proved due to significant and negative Spearman's coefficient, reflecting the moderate strength of relationship between them under energy range from 0.581 to 0.585.
So, below are the main inferences of our work, confirming in detail the abovementioned conclusion:

1.
In some cases, the single 2 H/ 1 H replacement result in positive value of OS occurrence frequency in the IFNA17 throughout E H cr diapason from 0.586 × 10 −22 N·m to 0.588 × 10 −22 N·m opposing to the OS occurrence frequency, when all hydrogen bonds in DNA nucleotides are 1 H, which is always equaled to 0.0 after E H cr exceed 0.585 × 10 −22 N·m. This underlines that at least 22.7% of the total number 2 H-substituted nucleobases can reduce molecular interaction in the studied gene and increase the hydrogen bond dissociation, foremost in its part from 328 to 653 nitrogenous bases; 2.
The counted occurrences of the OSs were much higher when the single 2 H/ 1 H replacement had taken place at nucleobases of the middle part IFNA17 (from 328 to 653 nucleotides) compared to its other parts, and this had a strong prevalence rate in the G-C pairs;

3.
The lowest rate of the OS occurrence was under the single deuterium substitution at the nitrogenous bases in the gene III part (from 654 to 980 nucleotides), which was also too rich in A-T pairs (72.2%) compared to the other parts of IFNA17, so that the calculated n max was equal to 0.0 for all of the studied critical energies (from 0.581 × 10 −22 N·m to 0.589 × 10 −22 N·m); 4.
Sum of n max was less significant when the single 2 H/ 1 H replacement occurred at the A-T nucleobase pairs compared to the G-C ones in the I and II parts of IFNA17, and the relationship between total n max and n max at G-C in these parts for each E H cr throughout the energy diapason from 0.581 × 10 −22 N·m to 0.588 × 10 −22 N·m was strong (r Spearman = 0.994) that proves the obvious dominant role of 2 H-substitutied G-C bases in the generation n max compared to the A-T; 5.
Earliest CSNBs (n = 3) arose under E H cr equal to 0.581 × 10 −22 N·m when IFNA17 had had the single 2 H-substituted cytosine or guanine nitrogenous bases in its I part (from 1 to 327 nucleotides). Moreover, throughout the range energy from 0.581 × 10 −22 N·m to 0.583 × 10 −22 N·m, the single 2 H/ 1 H replacement, leading to the CSNBs, prevailed in the I part, especially for its G-C pairs making up at least 67% of the total CSNB quantity in the whole gene. So, for IFNA17 throughout the range of E H cr from 0.581 × 10 −22 N·m to 0.585 × 10 −22 N·m the Spearman correlation coefficient between A-T/G-C ratio in the each gene part and CSNB numbers was significant and negative (r Spearman = −0.547); 6.
The highest acceleration of CSNB occurrence was observed when the single 2 H/ 1 H replacement took place at nucleotides of the III part of IFNA17 under E H cr from 0.585 × 10 −22 N·m to 0.586 × 10 −22 N·m, and throughout it they very abruptly reached the value of 100% of the nucleobases in this gene part. It indicates the obvious and higher vulnerability of IFNA17 due to the single 2 H-substitution at nucleobases from 654 to 980 compared with other gene parts exposed to studied critical energies, which increased the risk of permanent disorders of converting genetic information to mRNA messenger; 7.
All of the above-mentioned underline clearly the significant difference in the responsiveness of each IFNA17 parts under range of critical energies because of the single 2 H/ 1 H replacement in their nucleobases and with its strong dependence on A-T/G-C ratio with the prevalent contribution of the last pair, which leads to an increase due to the 2 H-substitution into its nitrogenous bases both n max and CSNB, especially under E H cr diapason from 0.581 × 10 −22 N·m to 0.584 × 10 −22 N·m; 8.
The single 2 H-substitution at G-C pairs not only had the most influence on the initial arise both n max and n CSNB in the whole gene throughout the critical energy diapason from 0.581 × 10 −22 N·m to 0.584 × 10 −22 N·m but also made possible the existence of the last OS occurrence under E H cr equals to 0.588 × 10 −22 N·at least in six cases that proves the leading effect of the isotopic 2 H/ 1 H modifications at G-C compared to A-T on the molecular dynamics of IFNA17; 9.
In addition, in the study was presented a modified algorithm allowing for accounting for nucleobases with the single 2 H/ 1 H replacement, which leads to occurrence of both the highest rate of Oss and CSNBs. Also, it showed the developed approach, decreasing significantly the false positive results compared to non-modified BJ-algorithm [29] due to the differentiated counting of the total sum of CSNB occurrence in the gene with relevance to the critical energy in the highest diapason.

Mathematical Model
The open states of the DNA molecule and its dynamics are well described by a mechanical model, which is two chains of disks that are interconnected by transverse springs; this system is described by the following Newton equations: here: ϕ i j (t)-is the angular deflection of the i-th nitrogenous base of the j-th chain counted counterclockwise at time t; I i j -is the rotational inertia of the i-th nitrogenous base of the j-th chain; R i j -is the distance between the center of inertia of the i-th nitrogenous 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 nitrogenous base pairs; F i j (t)-external influence on the i-th nitrogenous base of the j-th chain at a time t; n-is the number of nitrogenous base pairs in the system.
The magnitude of the external impact is assumed to be equal to F i j (t) = −β i j dϕ i j dt (t) + +F 0 cos ωt, where the term is −β i j dϕ i j dt (t) simulates the effects of dissipation caused by interaction with the water surrounding the DNA molecule, the term F 0 cos ωt is an external periodic effect.
In Equations (1)-(6), the first term to the right of the equality sign describes the force acting on the i-th nitrogenous base from the sugar-phosphate filament, the second term-the force from the complementary nitrogenous base, the third term-external impact.
Thus, Equations (1)-(6) allow us to model the hydrogen bond in the i-th pair (δ i = 1, k i 12 = k H,i 12 ), deuterium (δ i = 1, k i 12 = k D,i 12 ) and the break of this connection (δ i = 0). We will assume that a break in base pairs occurs if the potential binding energy in these pairs exceeds a certain critical value E H cr for a hydrogen bond and E D cr for a deuterium bond, if the potential energy in a pair with a broken bond is less than the critical value, then the bond is restored.
To Equations (1)-(6) we add the initial conditions: For the sake of certainty, we assume that at t = 0 the system is in equilibrium, that is, in the initial conditions (7) and (8) Problems (1)-(8) is a Cauchy problem for a system of 2n ordinary differential equations; in this paper, all studies were carried out on the basis of a numerical solution of this system.
The The value of the coefficient will be chosen because the deuterium bond is 5% stronger than the hydrogen bond. The order of the critical energy E is consistent with the experimental data from the work.
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 the i-th nitrogenous 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)- (8). To do this, we will create a set of points tj = jƮ, j = 1, , Ʈ = T/m in the segment [0, T]. Calculate at t = tj the ratio qj of the number of base pairs with a broken bond to the total number of base pairs n, then the value of Pk is equal to the arithmetic mean value over the points tj of these ratio:

Modification of Basov-Jimack Algorithm
The affiliation of each nucleotide to I part, II part or III part of IFNA17 was determined due to its sequence (serial) number. Nucleobases with , and the higher P were selected for the maximum range and their sum was , and nucleobases, which had CSNB (closed state of nitrogenous bases) ( = 0), were selected for the minimum range, and their number was designated as CSNB . According to the below-presented modified Basov-Jimack algorithm (BJ-algorithm [30]), which was changed for the highest critical energy range, every was arranged from CSNB to and their numbers were calculated in each gene part: where is number of i, which were included in the range "Maximum"; CSNB: CSNB = 0; E is critical energy in the diapason, which has 1 or more CSNB: 0 < CSNB / ≤ 1, in this case from 0.581 × 10 N·m to 0.589 × 10 N·m; is number of nitrogenous base pairs in x part of gene: x can be I, II or III part of IFNA17; CSNB is = 0.4 × 10 12 s −1 . We assume that E D cr = k D ·E H cr , k D,i 12 = k D ·k H,i 12 if one of the hydrogen bonds in the i-th base pair replaced with deuterium, k D = 1.05. The value of the k D coefficient will be chosen because the deuterium bond is 5% stronger than the hydrogen bond. The order of the critical energy E H cr is consistent with the experimental data from the work. 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 the i-th nitrogenous base pair one <any> hydrogen bond is replaced by deuterium.
The probabilities P 0 и P i , i = 1, n, will be sought on the basis of a numerical solution of the problem (1)- (8). To do this, we will create a set of points t j = j The value of the coefficient will be chosen because the deuterium bond is 5% stronger than the hydrogen bond. The order of the critical energy E is consistent with the experimental data from the work.
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 the i-th nitrogenous 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)- (8). To do this, we will create a set of points tj = jƮ, j = 1, , Ʈ = T/m in the segment [0, T]. Calculate at t = tj the ratio qj of the number of base pairs with a broken bond to the total number of base pairs n, then the value of Pk is equal to the arithmetic mean value over the points tj of these ratio:

Modification of Basov-Jimack Algorithm
The affiliation of each nucleotide to I part, II part or III part of IFNA17 was determined due to its sequence (serial) number. Nucleobases with , and the higher P were selected for the maximum range and their sum was , and nucleobases, which had CSNB (closed state of nitrogenous bases) ( = 0), were selected for the minimum range, and their number was designated as CSNB . According to the below-presented modified Basov-Jimack algorithm (BJ-algorithm [30]), which was changed for the highest critical energy range, every was arranged from CSNB to and their numbers were calculated in each gene part: where is number of i, which were included in the range "Maximum"; CSNB: Problems (1)-(8) is a Cauchy problem for a system of 2 ordinary differential equations; in this paper, all studies were carried out on the basis of a numerical solution of this system.
The study of the effect of 2 H/ 1 H exchange on the formation and dynamics of open states will be carried out using the example of the gene encoding interferon alpha 17. For this gene, = 980, the values of the coefficients of Equations (1)-(6) are taken from [31]), the values of = 0.526 × 10 , ɷ = 0.4 × 10 . We assume that E = • E , , = • , if one of the hydrogen bonds in the ith base pair replaced with deuterium, = 1.05. The value of the coefficient will be chosen because the deuterium bond is 5% stronger than the hydrogen bond. The order of the critical energy E is consistent with the experimental data from the work. 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 the i-th nitrogenous 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)- (8). To do this, we will create a set of points tj = jƮ, j = 1, , Ʈ = T/m in the segment [0, T]. Calculate at t = tj the ratio qj of the number of base pairs with a broken bond to the total number of base pairs n, then the value of Pk is equal to the arithmetic mean value over the points tj of these ratio:

Modification of Basov-Jimack Algorithm
The affiliation of each nucleotide to I part, II part or III part of IFNA17 was determined due to its sequence (serial) number. Nucleobases with , and the higher P were selected for the maximum range and their sum was , and nucleobases, which had CSNB (closed state of nitrogenous bases) ( = 0), were selected for the minimum range, and their number was designated as CSNB . According to the below-presented modified Basov-Jimack algorithm (BJ-algorithm [30]), which was changed for the highest critical energy range, every was arranged from CSNB to and their numbers were calculated in each gene part: (2) i ϵ range "Minimum" (CSNB): if 0 < CSNB / ≤ 1: CSNB = CSNB (E ); where is number of i, which were included in the range "Maximum"; CSNB: CSNB = 0; E is critical energy in the diapason, which has 1 or more CSNB: 0 < CSNB / ≤ 1, in this case from 0.581 × 10 N·m to 0.589 × 10 N·m; is number of nitrogenous base pairs in x part of gene: x can be I, II or III part of IFNA17; CSNB is number of closed states of nitrogenous bases in the gene part x under specific E ; is = T/m in the segment [0, T]. 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:

Modification of Basov-Jimack Algorithm
The affiliation of each nucleotide to I part, II part or III part of IFNA17 was determined due to its sequence (serial) number. Nucleobases with P i max , and the higher P i were selected for the maximum range and their sum was n max , and nucleobases, which had CSNB (closed state of nitrogenous bases) (P i = 0), were selected for the minimum range, and their number was designated as n CSNB . According to the below-presented modified Basov-Jimack algorithm (BJ-algorithm [30]), which was changed for the highest critical energy range, every P i was arranged from CSNB to P i max and their numbers were calculated in each gene part: if P i max − 1 10 P i max − P i min ≥ P 0 + 1 2 (P i max − P 0 ) and P i max > P 0 ≥ P i min ≥ 0: if P i max − 1 10 P i max − P i min < P 0 + 1 2 (P i max − P 0 ) and P i max > P 0 ≥ P i min ≥ 0: P i ≥ P i max − 1 4 (P i max − P 0 ) ⇒ n max = ∑ nP i ; (2) i range "Minimum" (CSNB): if 0 < n g CSNB /n g ≤ 1 : n CSNB = n x CSNB E H cr ; where n max is number of i, which were included in the range "Maximum"; CSNB: P CSNB = 0; E H cr is critical energy in the diapason, which has 1 or more CSNB: 0 < n g CSNB /n g ≤ 1, in this case from 0.581 × 10 −22 N·m to 0.589 × 10 −22 N·m; n x is number of nitrogenous base pairs in x part of gene: x can be I, II or III part of IFNA17; n x CSNB is number of closed states of nitrogenous bases in the gene part x under specific E H cr ; n g is number of nitrogenous base pairs in the whole gene.

Statistics
In addition, some statistical methods were used to exact the significance of differences among n CSNB and n max of the three gene parts. Yates corrected chi-squared test (χ 2 Yates ) was applied for a 2 × 2 contingency table (where degrees of freedom (ν) = (2 − 1)·(2 − 1) = 1). As a short-cut, for a 2 × 2 table with the following entries (Table 4): Table 4. 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 CSNB 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 CSNB of the range "Minimum", or n max of the range "Maximum" in the each compared part of the gene in the determined range of E H cr ; b and d equals the total number of nucleotide pairs of the gene part minus n CSNB of the range "Minimum", or n max of the range "Maximum" in the each compared part of the gene in the determined range of E H cr ; N A = a + b; N B = c + d; N S = a + c; N F = b + d; N = N A + N B + N S + N F .
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 CSNB of the range "Minimum", or n max of the range "Maximum" in the determined range of E H cr for two gene parts are mutually independent. The assess of the relationship between two variables was measured by Spearman's rank correlation coefficient (r Spearman ).

Conclusions
Thus, in the study, it was found that under the highest critical energy range, the very significant inequality of the frequencies of OS occurrence in the full gene were due to the single 2 H-substitution at the nucleobase of its different regions, which were almost equal to each other (gene I, II, III parts) by the number of nucleotides. In almost the entire highest critical energy range after the single 2 H/ 1 H replacement at the nitrogenous base of IFNA17, the P i values both in the range "Maximum" and "Minimum" were much different compared to the OS occurrence frequency under condition when all hydrogen bonds in DNA are 1 H (P 0 , Figure 2). Foremost, the A-T/G-C ratio had a strong influence on the frequency of CSNB in the E H cr diapason from 0.581 × 10 −22 N·m to 0.588 × 10 −22 N·m that can be not only the reason of the function impairment of IFNA17 but can also demonstrate the nucleic acid adaptation mechanisms to critical conditions and matter to the evolution process of living system [57]. For instance, under the highest critical energy range, gene function impairment can occur because of the increase in n CSNB , leading to the sharp slowdown of DNA bubble generation, which has more than four nucleobases, participating in forming of the preexisting denaturation state that is required by specific DNA-binding proteins [58,59]. In turn, via the single 2 H/ 1 H replacement at certain bases, the n max arising can be the mechanism of preserving the vital functioning of IFNA17 under critical conditions. The modified algorithm and new approach allow for calculating the n max and n CSNB throughout the highest critical energy diapason, in opposition to the non-modified Basov-Jimack algorithm, which is appropriate for counting the base-pair number in the range "Maximum" and "Minimum" under only the critical energies less than the highest range.
In addition, it should be noted, that the developed mathematical modeling showed its appropriateness for counting of the OS occurrence under the whole highest critical energy range. However, our study had a limitation, which was following: all results had been obtained in the frameworks of the mathematical model based on Newton's equations and are representative of the Cauchy problem for the system of 2n ordinary differential equations, which did not entirely take into account the all contribution of DNA-protein interactions [62,63]. Nevertheless, we do not contemplate this limitation as an insurmountable obstacle and consider the modified algorithm and new approach as appropriate methods that can be generalized and used for more other complex models of DNA dynamics.