Prof Piotr Cysewski

Ph.D.; Piotr Cysewski
email: piotr.cysewski@cm.umk.pl

PhD Beata Szefler

Environment influence on the aromatic character of nucleobases and amino acids

Keywords: aromaticity, nucleobases, amino acids, HOMA, NICS Aromaticity  Nucleobases  Amino acids  HOMA  NICS

Geometric (HOMA) and magnetic (NICS) indices of aromaticity were estimated for amino acids and nucleobases aromatic rings. The Cartesian coordinates were taken directly either from the PDB files deposited in public databases obtained with the finest resolution (≤1.5Å) or from structures resulting from full gradient geometry optimisation within hybrid QM/MM approach. Significant environmental effects imposing alterations of HOMA values were noticed for all analysed aromatic rings. Besides, even extra fine resolution (≤1.0Å) is not sufficient for direct estimation of HOMA values based on Cartesian coordinates provided by PDB files. The values of mean bond errors seem to be much higher than 0.05Å as it is often reported for PDB files. The quantum chemistry geometry optimization is strongly advised to be used and even simple QM/MM model comprising only the aromatic substructure within QM region and the rest of biomolecule treated classically within MM framework proved to be a promising way for description of aromaticity inside native environments. According to the presented results there are observed three consequences of of the interaction with the environment , inducing changes in structural and magnetic indices of aromaticity. First of all usually broad ranges of HOMA or NICS values are obtained for different conformations of nearest neighbourhood. Next, these values and their means can significantly differ from those characterising isolated monomers. The most significant increase of aromaticities is expected for six-membered rings of guanine, thymine and cytosine. The same trend is also noticed for all amino acids inside proteins but this effect is much smaller and reaches the highest value for five-membered ring of tryptophan. Explicit water solutions impose similar changes on HOMA and NICS distributions. Thus, environment effects of protein, DNA and even explicit water molecules are non-negligible sources of aromaticity changes appearing in rings of nucleobases and aromatic amino acids residues.

Response to Reviewers: 1.

geometries in all considered cases? I my opinion, the authors should clarify the following points (or at least adequate references should be given): 1) charge transfer effect: especially for the H-bonded ringsshould not be a QM region enlarged in these cases? 2) influence of the solvent: especially for the solvent exposed rings-should not be optimizations performed in the implicit solvent? 3) how correct pKa values of the ionizable groups in proteins were determined? 4)influence of the metal ions. Experimental structures do not offer positions especially of the loosely bound metal ions. How this problem was solved? …" I thank Referee#3 for pointing these warnings. In fact all these comments are related to methodology applied in this paper. First of all I apologize for not precise description of applied approach and most of posed questions come from this source. Thus, I have changed text describing method of QM/MM computations. Besides, since the paper was submitted for publication half year ago much more computations were performed in the meantime. For example explicit water solutions were studied using QM/MM method and collected data were provided in new figure in the revised version of this manuscript. This figure collects all obtained HOMA distributions in explicit water solutions and bioenvironment (proteins in case of amino acids and B-DNA for nucleobases). Thus, comments 1 and 2 are answered directly. I just emphasize that all HIGH regions in QM/MM studies comprised analyzed monomer along with all atoms within 3Å vicinity from any atom of the monomer. For explicit water solutions the inclusion of at least first hydration layer within QM region is important as it we previously demonstrated for fenol-fluoride system [Cysewski et al., , New J. Chem. 2009, 33, 831]. In case of proteins or DNA, however only monomers were included in QM region. The nucleobases interaction via hydrogen bonds, inter-or intra-strand stacking would require significant increase in computation resources needed. For example in case of nucleobases the QM region would increase 6fold after inclusion of only closest contacts. However, such system are too large for computations and crude model was applied instead. This is of course the source of error but I believe it does not seriously affect purpose of this paper. There are two aims of our work. First of all there is offered critically analysis of the direct applicability of Cartesian coordinates provided by high quality experimental data collected in PDB files. The second goal is to provide evidences that the environment can significantly affect aromatic ring flexibility quantified by some aromaticity indices. This paper is not devoted to detailed analysis of particular neighborhood (as B-DNA or protein sequence, context or conformation). Just statistically significant probe of structures and non-zero (non-trivial) distributions are obtained, which are worth of its own. Of course more detailed studies are to be done and this will be the subject of our further investigations. In this context comments 3 and 4, although valuable and generally very important, are less related to submitted paper.

Introduction

Environment affects of physical and chemical nature can introduce significant structural alterations to molecules directly exposed to such forces. These changes are of particular interest for molecules important from the bio-chemical perspective. Nucleic acid bases and amino acids, as important building blocks of polynucleotide and polypeptide chains, are compounds of particular significance. All nucleobases and four amino acid residues posses aromatic rings, which can be affected by environment heterogeneity. These influences can be as different as explicit solvent molecules, dense packing in solids, particular molecular surroundings in polymeric neighborhoods and variety of organic or inorganic ligands directly interacting with solute. Such straightforward contacts between molecules in closest proximity, chemical modification or indirect influence via intermolecular interaction can impose sizable alterations of physico-chemical properties including aromaticity. This term is associated with a cyclic π-electron delocalization imposing an increase of stability measured by aromatic stabilization energy [START_REF] Minkin | Aromaticity and Antiaromaticity Electronic and Structural Aspects[END_REF][START_REF] Cyrański | [END_REF][3], small bond length alterations within ring skeleton [4][5][6], π-electron ring current formation after exposition to external magnetic field [7][START_REF] Lazzeretti | [END_REF] and typically reactivity toward substitution rather than addition reaction [9,[START_REF] Smith | Advanced Organic Chemistry[END_REF]. This proves the fact that aromaticity is a multidimensional phenomenon [START_REF] Krygowski | [END_REF] and justifies utilization of different measures of aromatic character. Most of aromaticity indices are not independent since in many reports there were documented correlations between different aromaticity criterions [12][13][14]. Despite the fact that the mutual relationships between aromaticity indices strongly depend on the selection of molecules in the probe [12,15,16] the correspondence between different indices can be significant. For example in the case of porphyrins [13] the magnetic index changes monotonically with the geometry-based index. Similarly, rings in benzenoid hydrocarbons exhibit linear correlation between HOMA and NICS [14]. Besides, energies of ring stabilization in benzenoid hydrocarbons correlate with NICS and HOMA values [17,18]. Moreover, magnetic index of aromaticity correlates best with geometry based indices for fivemember heterocyclic compounds [12]. Although not all aromaticity indices can be applied for heterocycles as for example A j , index formulated by Julg and Francois [19,[START_REF] Julg | Aromaticity, Pseudo-Aromaticity, Anti-Aromaticity[END_REF] or bond alternation coefficient (BAC) formulated by Krygowski and co-workers [START_REF] Krygowski | [END_REF] but there are many indices effectively used for this particular class of compounds. Especially HOMA index proved to be universal tool offering not only the excellent agreement with chemical experience and intuition [22] but also it can describe electron delocalization of fragment or the π-electron system as the whole. There is additional and very important advantage of using HOMA as the criterion of aromaticity, namely the potential applicability to experimental bond lengths resulting from variety of diffraction measurements. These data are deposited in public databases providing molecular geometries of organic molecules [23] (Cambridge Structural Database), amino acids [24] (Protein Bank Database) or nucleobases [25] (Nucleic Acid Database). The usefulness of this kind of applications was proven to be effective in numerous studies [26][27][START_REF] Domenicano | Accurate Molecular Structures: Their Determination and Importance[END_REF][START_REF] Colapietro | [END_REF][30][31].

Intermolecular interactions are important factor affecting aromaticity. For example most heterocyclic and some carbocyclic compounds increase their aromaticities with an increase of medium polarity [32]. Also differences of net of hydrogen bonding may impose seriously aromaticity alterations, as for example it takes place in the case of p-nitrosophenolate anion in crystalline state [33]. Furthermore, aromaticity of nucleobases involved in Watson-Crick pairs also can be significantly changed with respect of the isolated monomers [34]. The systematic study of the hydrogen bonding on aromaticity of phenol and p-nitrophenol were undertaken by Krygowski et al. [35] demonstrating monotonous relation between rings aromaticities and the H-bond strength. The perturbation of aromaticities imposed on π-electron system by alkali metal cations [36][37][38][39][40][41], alkaline earth cations [42][43][44][45] as well as transition metal cations interactions with benzene [46][47][48] were also intensively studied. Furthermore, aromaticities of nucleobases and amino acids were described based on different aromaticity criterions [49][50][51][52][53][54][55][56][57][58][59][60][61][62].

The aim of this study is to analyze the environmental influence on the aromatic character of four amino acid residuess (phenyloalanine (PHE), tyrosine (TYR), tryptophan (TRP) and histidine(HIS)) and DNA bases (adenine (ADE), guanine (GUA), thymine (THY) and cytosine (CYT)) by the aid of geometric (HOMA) and magnetic (NICS) indices. This offers evidences that the environment can significantly affect the aromatic ring flexibility, as quantified by these aromaticity indices. Besides, the direct utilization of molecular geometries taken from public databases is discussed and compared to QM/MM optimized structures of aromatic rings immersed in native environment. The separation of influences on aromaticit y of actual protein or DNA surroundings from artifacts related to experimental inaccuracies is another important aspect considered in this study. To our best knowledge this is the first study demonstrating sensitivity of HOMA index to experimental errors for considered aromatic compounds.

Methods

The Cartesian coordinates of amino acids were extracted from Brookhaven Protein Database [24,63] after the selection of only those structures which met quality standard of the finest resolution (≤1.5 Å). Similarly, the coordinates of nucleobases were retrieved from Nucleic Acid Database also using the above criterion of accuracy. The probe of aromatic rings is quite extensive and the statistical details of data distributions are provided in the supplementary materials (Tables S1 andS2 for two classes of the analysed compounds). Two criterions of aromaticity were applied, namely HOMA [64][65][66] and NICS [67] indices. These two criterions were selected due to their wide usage [13,[67][68][69], simple meaning [65,66] and direct applicability to both experimental [26][27][START_REF] Domenicano | Accurate Molecular Structures: Their Determination and Importance[END_REF][START_REF] Colapietro | [END_REF][30][31]70] and theoretical geometries [35,[71][72][73][74][75][76].

Three types of calculations were performed in this project to estimate the aromatic properties of analysed compounds. In the first step, the estimation of HOMA sensitivity on the experimental errors of bond lengths was performed. In the first step the model analysis was performed for estimation of HOMA sensitivity on the experimental errors of bond lengths.

Then HOMA and NICS(1) values were estimated directly using experimental geometries.

Finally, these two criterions of aromaticity were calculated for geometries resulting from quantum chemistry optimisations within QM/MM hybrid approach [65,66]. The model analysis was performed for generation of HOMA distributions related to geometry distortions induced by experimental inaccuracies. The original [65,66] formula defining HOMA was modified by introduction of the following bond length changes:

       i j i opt j R R n HOMA 2 ) ( 1 1 ) (   
where ω j stands for the error of i-th bond length (in case of ω j =0 classical definition of HOMA is recovered), R opt corresponds to reference values [66] (1.388 for C-C bond and 1.334 for C-N bond), R i represents the optimised bond length for isolated monomer in the gas phase, α is normalisation constant [66] equal to 257.7 and 93.52 for C-C and C-N bonds, respectively. The ω j values were taken in the range of <-0.1 Å,+0.1 Å> with 0.005 Å step and were used for the alterations of bonds lengths. All possible combinations of bonds deformations were taken into account by successive changes of one, two, three, etc, up to all bonds constituting the ring. For example in case of six membered ring there are 704 unique combinations for any value of bond length errors. The five membered ring may be affected on up to 19008. On such a way series of HOMA values were generated corresponding to all possible distortions of the optimised structures by increasing or decreasing the bond lengths in the ring skeleton Then structures were grouped according to the distribution of the average bond length errors, as shown by the smoothed histograms in Fig. 1, 5 and 6. In the second part of our study the full gradient optimizations were performed for representative sets of aromatic rings. The hybrid quantum chemistry/molecular mechanics approach was used for systems containing particular aromatic amino acid residue or nucleobase with initial geometries taken directly from PDB files. The geometry of particular residue (PHE,TYR, HIS, TRP or any nucleoside) was optimized on B3LYP/6-311+G** level of theory. The rest of the system was treated classically using AMBER force-field [79,80] embedded according to ONIOM [77,78] methodology. Besides, explicit water solutions were modeled by immersing of nucleobases or aromatic amino acid residues within the drop of water of 20 Å diameter. The QM part in these QM/MM studies comprised analyzed monomer along with all water molecules within 3 Å vicinity from any atom of the monomer. Inclusion of at least first hydration layer within QM region is important as it was previously demonstrated for fenolfluoride system [75]. In case of proteins or DNA, however only monomers were included in the QM region. The nucleobases interaction via hydrogen bonds, inter-or intra-strand stacking would require significant increase in computation resources needed. For example in case of nucleobases the QM region would increase 6-fold after inclusion of only closest contacts. This poses serious technical problems due to rapid increase of computational resources needed. Although such simplified model can be the source of error due to ignoring important influences (mainly dispersions) but is proposed here as a compromise between the accuracy and the cost of computations. Some tests were performed for amino acid molecules inside protein and inclusion of the nearest neighboring amino acids into QM part only slightly affected HOMA values (with increasing HOMA by <0.01 unit). The extension of QM region led to significant increase of computations costs and difficulties in convergence of SCF procedure and minimization protocol. Thus, the proposed model properly serves to the purposes of the paper. All quantum chemistry and hybrid QM/MM calculations were performed using Gaussian03 package [81].

Results and discussion

Definition of the HOMA index accounts for its significant sensitivity on the bond length deviations from the reference values. The intention of using the square-like relation was the quantification of aromaticity changes induced by even minute alterations from ideal structure of benzene. However, if coordinates are directly taken from X-ray diffraction patterns of crystals the inaccuracies coming from the experimental and interpretation procedures can also affect geometric criterion of aromaticity. These changes naturally are not necessarily related to the alteration of aromatic character but can just represent artefacts of inaccurate geometries.

Only in the cases were observed variations of HOMA are higher than those related to individual bond length errors, the environment influence on the aromatic character can be considered as significant. In order to visualise the influence of experimental inaccuracies on HOMA values the series of model calculations were performed as it was described in the methodology section. Since different compounds herein analysed have different aromatic nature and different rings size, the analysis was performed separately for each of the class of aromatic rings.

Experimental errors affecting HOMA values of amino acids

In Fig. 1 Consequently, there is fairly linear relationship between skewness and kurtosis estimated for different bond length inaccuracies as it is shown in Fig. 2. Since for normal distributions kurtosis and skewness are equal to zero these quantities can be used as an indicator of HOMA alterations observed for experimental geometries. It is evident that the most sensitive to the geometry inaccuracies are aromatic rings of phenyloalanine and tyrosine since comparable bonds distortions lead to the highest values of kurtosis and skewness. Tryptophan as compound comprising two conjugated rings is much more resistant to experimental errors of ring deformations. Also aromaticity of five-membered rings are slightly less prone to bond length errors compared to six-membered rings. In Fig. 1 there are also presented plots of smoothed histograms obtained for aromatic rings of amino acid residues, which are directly exposed to influences of protein environment. The detailed statistical characteristics is also provided in Table S1 (see supplementary materials). All PDB files used for Cartesian coordinates extraction belong to the best available structures obtained with resolutions better than 1.5 Å. Despite of this fact, surprisingly high values of skewness and kurtosis are observed for HOMA distributions corresponding to experimental geometries. As it is documented in Table S1 and Fig. 2a all values of skewness are negative and kurtosis can reach values of several units. Thus, either individual bond length suffer from more than 0.1 Å or environmental effects are strong and non-negligible. It is interesting to separate these two effects in order to get deeper inside into observed heterogeneity of HOMA values. First of all, we have found that about 2% of experimental structures were obviously incorrect. Among several thousands of rings that were analysed, some of them were incomplete for not having all atoms in the PDB files. In other cases geometries were incorrect since their use lead to non-physical values of HOMA (for example negative values for PHE). Thus, simple estimation of HOMA index based on experimental geometries can be alternative, direct and effective way of checking against eventual errors of the Cartesian coordinates obtained via Xray diffraction protocols. It is claimed that for files obtained with fine resolution (≤1.5 Å) the estimated overall coordinate error is usually less than 0.05 Å. However, from the presented plots in Fig. 1 it is obvious that at least for aromatic rings the actual distortions must be much higher. This is also addressed in Fig. 3, where the correspondence between the HOMA values estimated using Cartesian coordinates taken directly from PDB files (denoted hereafter HOMA PDB ) and those obtained from the QM/MM optimisation (HOMA QM/MM ), is presented.

It is evident that the broad range of HOMA PDB values is significantly narrowed after QM/MM geometry optimisation. In case of tyrosine the QM/MM optimisation did not introduce any changes for ring aromaticities since all HOMA PDB and HOMA QM/MM values differ by less than 0.012 units. The six-membered ring of tryptophan proved to be slightly more sensitive to optimization with corresponding changes from 0.068 to 0.023 units. For all other amino acid residues the geometry was significantly changed by QM/MM optimizations. For example for phenylalanine the discrepancies in HOMA values were reduced from 0.170 to 0.017 units of HOMA. The most significant impact of the QM/MM optimisation is observed for five-membered rings of histidine (reduction from 0.576 to 0.079) and tryptophan (decrease from 0.615 to 0.072). This suggest that most of HOMA heterogeneity coming from PDB files should not be directly related to protein environment but rather to experimental inaccuracies. This is also supported by the fact that the localisation of maxima in Fig. 1. indicating the most frequently occurring values of HOMA are almost the same for PDB structures and those coming from model analysis. This is also supported by the fact that the localisation of maxima in Fig. 1 are almost the same for PDB structures and those coming from model analysis. This conclusion is not surprising in the light of suggesting that the πelectron delocalisation in aryl ring is strongly resistant to chemical modifications of substituents. However, it is known [35,[72][73][74][75][76] that environment can significantly affect aromaticity as for example in the case of fluoride approaching phenol in the gas phase leading to significant reduction of the aromaticity of the phenyl ring [35,75]. Although this effect is significantly screened [75] by the presence of discrete water molecules it is still nonnegligible, particularly emphasized by the broad spectrum of HOMA values obtained for different conformations of the first hydration layer around the phenol solute. Data presented in Fig. 3 suggest that protein environment only slightly affects amino acid aromaticities. However, systematic trends are observed, the HOMA values of the most of QM/MM structures being higher in comparison to those of the isolated monomers in the gas phase. The five-membered rings of histidine and tryptophan are characterised by the highest diversity of HOMA values obtained after QM/MM optimisation. Furthermore, in case of histidine, interesting differences in the shape of HOMA distribution plots corresponding to fine (≤1.5 Å) and extra fine (≤1.0 Å ) resolution (see Fig. 1) are observed. One obvious reason of this fact might be the difference in the size of both samples. However, for other amino acids this effect is almost insignificant. The second observation is that the set of highest resolution histidine residues is not quite representative since it comprises significant portion of histidineporphyrine complexes. Among 47 histidine residues found within PDB files corresponding to extra fine resolution (≤1.0 Å) there are 16 rings that form complexes with HEM moiety via iron cations. These interactions might affect aromaticity of histidine. In order to quantify the presence of iron-histidine interactions the model system was analysed as shown in Fig. 4. It contains four amino groups and two histidine molecules forming complex with one Fe +2 cation. For such a simplified model of HEM the separation distance of one histidine molecule (HIS(1)) was fixed at specific values and all the remaining coordinates were optimised on B3LYP/6-311+G** level both in the gas phase and in the water solutions using PCM model. Interestingly, the separation distance of both histidine molecules from iron cation are not independent and the increase of r 1 bond forces the decreasing of r 2 values, what suggests approaching HIS(2) closer to metal centre. The shorter Fe-HIS bond the smaller aromaticity of histidine and this relationship is almost linear (with R 2 =0.980). The polar environment modeled by continuum field of water significantly reduces the influence of iron cation on the aromaticity of histidine. Despite the fact that Fe-HIS(2) bond length are slightly shorter in water solution the HOMA values are significantly less affected by this kind of environment interactions. The polar solution screens down structural changes induced by the presence of iron cation in the model of HEM system. This suggests that the observed differences of HOMA distributions between extra fine and fine resolutions is only partly explained by this kind of interactions. Indeed the average value of bond length formed between Fe 2+ and the nitrogen atom of histidine ring in PDB files is equal to 1.974±0.028 Å. However, closer inspection of PDB files leads to another interesting observation, namely there is significant difference between aromaticity of histidine residues place inside the protein and in the outer regions. The most aromatic histidine rings occupy external regions of protein. These structures are exposed to interactions with solution environment and can be more easily protonated, interact more directly with the water molecules or cations that are present in the solution. Thus, these external histidine moieties are responsible for increase of aromaticity and consequently affecting the HOMA distributions. On the other hand the histidine-HEM complexes are simple forced to be inside the protein volume and are prone to local interactions, which usually reduce aromaticity of this five-membered ring. Both these effects cause the flattening of HOMA values distributions. In the case of much larger set of structures corresponding to fine resolution (≤1.5 Å) the populations of such external localisations is significantly smaller and even the reduction of resolution to ≤2.0 Å of PDB files does not alter the shape of HOMA distribution. In conclusion it is worth emphasizing that although it is possible to identify some environmental effects that impose alterations of HOMA of aromatic amino acid rings the majority of observed HOMA variations are to be related merely to experimental inaccuracies rather than to protein environment itself and even extra-fine resolution (≤1.0 Å) is not sufficient for direct assessment of aromaticities via the geometric index.

Experimental errors affecting HOMA values of nucleobases

Analogue investigations were performed to estimate the environmental effect on aromaticity of nucleobases. Again model analysis was followed by extraction of the Cartesian coordinates from PDB files [24,63] with extra fine (≤1.0 Å) and fine (≤1.5 Å) resolutions. As illustrated in Fig. 5 and6, the main consequence of bond length distortions is the broadening of HOMA values distribution, without significant affecting the positions of maxima of smoothed histograms. This observation is the same as mentioned above for amino acids residues and the linear relationships between skewness and kurtosis is also fulfilled for nucleobases (as it can be seen on Fig. 2b). However, the analysis of HOMA distributions coming from PDB files reveals new interesting feature. First of all the aromaticity of pyrimidine rings in terms of the HOMA index is much higher inside polynucleotide chain than for isolated molecule. The only exception form imidazole and purine rings of adenine. On the contrary according, to HOMA PDB values the whole molecule aromaticity of guanine is increased by 0.15 units if this nucleobase is present within polynucleotide chains. Interestingly, the analysis of different polymorphic forms as double stranded (ds)A-DNA, B-DNA, Z-DNA or even single stranded (ss) RNA interiors does not alter these conclusions. The differences in HOMA values distributions in these various environments are rather minor, although some trends are visible.

For example maxima on smoothed histograms are located systematically at higher values of HOMA for ss-RNA compared to ds-B-DNA irrespectively of the nucleobase. Of course thymine is not present in the case of RNA and instead the uracil moiety is analyzed in Fig. 5.

An even higher impact of Z-DNA environment is observed on aromaticities of cytosine populations. Interestingly, this effect is not observed in the case of guanine, for which all analysed forms of polynucleotide chains induce very similar variations on HOMA distributions irrespectively of the ring size. Furthermore, as it was demonstrated in Fig. 2b the skewness and kurtosis of HOMA PDB values distributions are much higher for purine than for pyrimidine rings. They are even higher than for distributions obtained from model analysis. This might suggest that these rings are very sensitive to bond lengths errors or the environment effect are significant. The analysis presented in Fig. 7 clarifies this observation.

After optimisation of nucleobases using QM/MM approach, the variability of HOMA QM/MM is significantly reduced if compared to HOMA PDB . The most spectacular change is clearly visible in case of thymine, for which aromaticities of few conformations are completely incorrect if PDB coordinates are used. After geometry optimization structures adopt right geometry univocally addressing strange values of HOMA PDB as coming from coordinates inaccuracies. Also cytosine is prone to experimental inaccuracies. After nucleoside geometries optimisations the smallest variations of HOMA QM/MM values are observed for adenine (0.024) and the highest for cytosine (0.106). Furthermore significant increase of rings aromaticities are observed inside DNA if compared to monomers in the gas phase. Only imidazole rings and purine rings of adenine seem to behave on the opposite manner. Thus, aromaticity can be significantly affected by different context and structure of polynucleotide chain leading to considerable alteration of HOMA values with respect of free monomers. The importance of the environmental effects of nucleic acid interiors on the aromaticities of nucleobases is more pronounce than for protein environment on amino acid residues. The observed quite significant increase of both HOMA PDB and HOMA QM/MM vales should not be related merely to errors imposed on coordinates by X-ray diffraction processing, although they cannot be neglected. Furthermore the polymorphic forms are important factors affecting nucleobases aromatic characters and the direct use of Cartesian coordinates significantly exaggerate this effect. Theses aspects deserve a more detailed analysis and they will be addressed in the subsequent study.

NICS(1) alterations imposed by environment

The final aspect analysed in this paper is the influence of the environment on aromaticity expressed by magnetic susceptibility index. The negative values of magnetic shielding estimated 1 Å above or below the ring centre (NICS(1) or NICS(-1)) is commonly used as a measure of the magnetic criterion of aromaticity. The geometries obtained after QM/MM optimisations were used to estimate the NICS values of both amino acid resides and nucleobases. Obtained values along with ones corresponding to isolated monomers (optimised in the gas phase) were presented in Fig. 8. Interestingly six-membered ring of TRP as the most aromatic in gas phase among all amino acid residues is also the most affected by protein environment. Heterogeneity of aromatic character is quite high and amounts 1about 8% for TRP (5) and 7% for histidine. The lowest influence of the environment is observed for rings of phenyloalanine, tyrosine and six-membered rings of tryptophan. In the latter case it reaches 0.04 of NICS units with 4% of fluctuations. For the most conformations comprising aromatic rings of amino acids, the reduction of NICS values with respect of the isolated monomers is noticed. The most significant influence of protein environment on NICS values is observed for five-membered rings of tryptophan and histidine. The mean values of ΔNICS estimated as difference with respect to the isolated monomers (NICS protein -NICS gas ) are equal to -1.0ppm and -0.5ppm for TRP (5) and HIS (5), respectively. On the contrary the averaged aromaticity of six-membered ring of tryptophan is almost not affected by protein environment (ΔNICS=0.0).

Phenylalanine and tyrosine are characterised by ΔNICS equal to -0.4 and -0.1 respectively. All these conclusions are in good accord with previously discussed observations related to the structural index of aromaticity. Thus, the protein environments modestly affect the aromaticity of amino acids, although for some conformations such influence is nonnegligible.

Similar results, characterising the B-DNA influence on the magnetic index of nucleobases aromaticity is presented in the part b of Fig. 8. In case of double stranded helix, two contexts corresponding to NICS(1) and NICS(-1) distributions are clearly identified. Here, the former is related to 5'-side and the latter to 3'-region of nucleobase. Since the point located above purine rings is above C-C bond the corresponding NICS values cannot be used as a measure of the whole molecule aromaticity and only pyrimidine and imidazole rings are analysed here.

According to this measure of aromaticity, the polynucleotide chain can significantly affect aromaticities of all rings since the span of obtained values ranges from the the smallest variability equal to 1.1ppm in case of 3'-imidazole ring of adenine up to the highest variety for 3'-thymine rings reaching 3.3ppm. Both analysed aromatic rings of adenine and imidazole ring of guanine are characterised by positive values of mean ΔNICS, what suggests a decrease of aromaticities with respect of the free monomers in the gas phase. For the rest of rings significant increase of aromaticity is expected. For example the mean ΔNICS values of 5'and 3'-pyrimidine rings of guanine are equal to 1.1ppm and 2.0ppm, respectively. Also the rest of 5'-(3'-) pyrimidine rings are characterised by increase of their aromatic character by amount of -0.5(-0.8)ppm and -1.4(1.0)ppm for cytosine and thymine, respectively.

Conclusions

Presented above discussion demonstrated that even extra fine resolution (≤1.0 Å) is not sufficient for direct estimation of HOMA values based on Cartesian coordinates provided by PDB files obtained via X-ray diffraction of protein and DNA crystals. The values of mean bond errors seem to be much higher than 0.05 Å as it is often reported for such files. This imposes serious limitations on direct use of Cartesian coordinates and even the best available PDB files cannon provide structures of acceptable accuracy. On the other hand the estimation of HOMA values is proposed as an effective and alternative way of verification of accuracy of experimental structures. Despite the fact that environment effects are significantly overestimated if direct coordinates are taken from the PDB files they are still are nonnegligible. However, quantum chemistry geometry optimisation is strongly advised to be used and even simple QM/MM model comprising only the aromatic molecule within QM region and the rest of biomolecule being treated classically within MM framework proved to be a reliable way for aromaticities description inside native environment. According to presented results there are observed three consequences of environment inducing changes of structural index of aromaticity. First of all usually broad range of HOMA values are obtained for different conformations. These values and their mean values can significantly differ from those characterising isolated monomer. Also the protein and DNA interiors can affect the values of magnetic index of aromaticity both by increasing the range of values related to different conformations and altering the mean ΔNICS values estimated with respect of isolated monomers in the gas phase. The most significant increase of aromaticity is expected for pyrimidine rings of guanine, thymine and cytosine. The same trend is also visible for all amino acids inside proteins but this effect is much smaller reaching the highest value for fivemembered ring of tryptophan. Thus, the environment effects are non-negligible and need further investigations to reveal the details of aromaticity changes imposed by biomolecule interiors Interestingly, also explicit water molecules form the significantly non-homogeneous environment. In

Figure captions

Fig. 1 Heterogeneity of the structural index of aromaticity (HOMA) calculated for four aromatic amino acids residues. Two sets of plots correspond to systematic analysis of mean bond length errors (gray lines) and to protein environment (black lines with open circles or triangles). Only high quality structures (resolution better than 1.5 Å) were taken from Protein Database [24]. For detailed statistics see Table SI.

Standard deviations of the analysed HOMA values distribution are given in brackets. The amount of mean bond length errors is indicated in the legends.

Symbol HOMA(x) denotes the aromaticity of the x-membered ring Heterogeneity of the magnetic measure of aromaticity (NICS(1)) of amino acids and nucleobases inside protein (PDB [24]) or DNA (NDB [25]) interiors, respectively. The lines represent reference levels of isolated monomers optimized in the gas phase

Fig. 9

Heterogeneity of the structural index of aromaticity (HOMA) corresponding to nucleobases or amino acids inside DNA (NDB [25]), protein (PDB [24]) interiors.

Besides, the gray symbols represent distributions of HOMA values estimated for rings geometries optimized in explicit water solutions using QM/MM approach.

The lines represent reference levels of isolated monomers optimized in the gas phase Fig.

1 6) GUA( 6) THY( 6) CYT( 6) ADE( 5) GUA( 5) ADE(9) GUA(9) Fig.

2

  distributions of HOMA values were presented for each of analysed aromatic ring belonging to the amino acid residue for three selected values of mean error of bond lengths. Significant influence of the inaccuracies on HOMA values is clearly visible irrespectively of the residue type. First of all there are considerable negative skew characterizing all distributions except only the most precise case when the bond length inaccuracies do not exceed 0.01 Å. Geometries of such accuracy lead only to small HOMA alterations what is expressed by very small values of standard deviations, skewness and kurtosis. The increase of the mean value of the bond length significantly affects HOMA distributions what is clearly visible by the presence of a strong negative skew. In all cases the highest values of mean error of bond lengths are accompanied by the broadest distributions of HOMA values.

  Fig. 9 collects values corresponding to QM/MM computations of biopolymer environments and water solutions. As it was discussed previously both B-DNA and protein interiors significantly affect aromaticity of the analysed compounds. What is also worthy noticing, the fluctuations of water molecules around aromatic rings impose flexibility of amino acid residues and nucleobases rings quantified by HOMA distributions. Although this paper is not devoted to detailed analysis of particular neighborhood (as B-DNA or protein sequence, context or conformation) the provided statistically significant probe of structures with non-zero (non-trivial) distributions of aromaticity indices offer serious evidence that aromaticity can be affected by environmental factors that do not induce chemical modifications of the analyzed aromatic compounds. Further investigations are indispensible for a more precise description of these contributions. W, Wong MW, Gonzalez C, Pople JA (2004) Gaussian 03, Rev C.02. Gaussian Inc, Wallingford CT 82. Krygowski TM, Stępień BT (2004) Pol J Chem 78:2213-2217 83. Krygowski TM, Stępień BT, Cyrański MK, Ejsmont K (2005) J Phys Org Chem 18:886-891

Fig. 2

 2 Fig. 2The relationship between skewness and kurtosis for HOMA distributions of amino acid (a) and nucleobases (b) aromatic rings affected by various bond length inaccuracies. Larger symbols represent skewness and kurtosis of experimental values estimated from Cartesian coordinates directly taken from PDB files[24,25]

Fig. 3 Fig. 4 Fig. 5 Fig. 6 5 ÅFig. 7 Fig. 8

 3456578 Fig. 3Correlation between HOMA values of aromatic amino acids estimated using Cartesian coordinates taken directly from PDB files (HOMA PDB ) and ones obtained after QM/MM optimisation (HOMA QM/MM ) at ONIOM(QM:B3LYP/6-311+G**//MM:AMBER) level. The vertical lines represent HOMA values estimated for geometries of isolated monomers optimized in the gas phase using B3LYP/6-311+G** method

  Figure 1 Click here to download line figure: Fig1.doc

  Figure 2 Click here to download line figure: Fig2.doc

  Figure 3 here to download line figure: Fig3.doc

Figure 4

 4 Figure 4 Click here to download line figure: Fig4.doc

Fig. 5

 5 Fig. 5

Figure 5

 5 Figure 5 Click here to download line figure: Fig5.doc

Fig. 6

 6 Fig. 6

  Figure 6 Click here to download line figure: Fig6.doc

  Figure 7 Click here to download line figure: Fig7.doc

  Figure 8 Click here to download line figure: Fig8.doc

Figure 9

 9 Figure 9 Click here to download line figure: Fig9.doc

Acknowledgments

The work was supported by the Computational grant no 39 of PCSS (Poznań, Poland). The allocation of computational facilities is greatly acknowledged.

Supplementary material available

 6) TYR( 6) TRP( 6) HIS( 5) TRP( 5) TRP(9) PDB files (a) 264 337 462 469 462 462 Residues (b) 801


Hui Kong
email: kong-h08@mails.tsinghua.edu.cn

Hehua Zhang

Xiaoyu Song

Ming Gu

Jiaguang Sun

Proving Computational Geometry Algorithms in TLA+2

Keywords: loop invariant, algorithm verification, theorem proving, TLA+2, geometry algorithm

Geometric algorithms are widely used in many scientific fields like computer vision, computer graphics. To guarantee the correctness of these algorithms, it's important to apply formal method to them. In this paper, we propose an approach to proving the correctness of geometric algorithms. The main contribution of the paper is that a set of proof decomposition rules is proposed which can help improve the automation of the proof of geometric algorithms. We choose TLA+2, a structural specification and proof language, as our experiment environment. The case study on a classical convex hull algorithm shows the usability of the method.

I. INTRODUCTION

Computational geometry is an important discipline in computer science. It mainly focuses on studying algorithms for solving geometric problems [START_REF] Preparata | Computational geometry: an introduction[END_REF]. The applications of geometric algorithms range over a variety of scientific fields such as computer graphics, computer vision, robotics, astrophysics. Therefore, guaranteeing the correctness of these algorithms is crucial, especially when they are used in safety-critical cases.

One common feature of computational geometry algorithms is their involvement in plenty of operations on large numbers of points, lines, or planes. Consequently, array is the most commonly-used data structure to store these geometric elements. This usually means that the loop invariants of these algorithms contain some predicates with quantifiers. Therefore, to prove the correctness of these algorithms, we need some approaches to proving the loop invariants of this form.

However, currently, there are no effective automatic methods for this purpose and it mainly depends on sophisticated researchers to do this work manually in interactive theorem provers like Coq etc [START_REF] Pichardie | Formalizing convex hull algorithms[END_REF]. This may work well for those illustrative simple examples, but for those complex realistic problems such as array-based geometric algorithms, the work will become error-prone and uncontrollable.

In this paper, we propose an approach to improving the automation of the proof of the loop invariants of geometric algorithms. Our interest lies in those whose loop invariants are predicates with quantifies. The basic idea is as follows:

1) structural decomposition. Given a geometric algorithm G and its candidate loop invariant I as well as the proof obligation G ∧ I ⇒ I , by transforming G and I into DNF(Disjunctive Normal Form) and I into CNF(Conjunctive Normal Form) respectively, the proof obligation is decomposed into a number of simpler proof obligations such that each of them is an implication whose precondition and consequence are a conjunction and a disjunction respectively. 2) quantifier decomposition. Based on the result of the previous decomposition, the consequence part of every sub-proof obligation should be a PNF(Prenex Normal Form) of the form

Q 1 x 1 ...Q n x n : x 1 ∈ D 1 ∧ ... ∧ x n ∈ D n ⇒ P (x 1 , ..., x n )
, where Q i is a universal or existential quantifier. According to the change to domains of quantified variables, each domain D i is divided into two parts: the unchanged part and the added part with respect to its corresponding domain D i before the loop. Based on this division, every sub-proof obligation is further decomposed into lots of much simpler cases which is often obvious to an interactive theorem prover or need a little human intervention.

As convex hulls are widely used in computational geometry, we choose a typical two dimension convex hull algorithm-Graham's Algorithm for our case study. We choose TLA+2 as the algorithm specification and proof language and TLA+ toolbox as the experiment environment.

The remainder of this paper is organized as follows. Section II briefly introduces the preliminaries of our method. Section III introduces the formal specification of convex hull and the Graham's algorithm in TLA+2. Section IV demonstrates on how to prove a loop invariant for a geometric algorithm. Section V explain briefly our experiment result. Section VI introduces related work. Finally, we conclude our work in Section VII.

II. PRELIMINARIES

TLA+2 is a formal specification and proof language based on TLA(Temporal Logic of Actions), a logic developed by Leslie Lamport, which combines temporal logic with a logic of actions [START_REF] Lamport | Specifying systems: The TLA+ language and tools for hardware and software engineers[END_REF].

TLA+2 is used to describe behaviors and properties of sequential or concurrent systems. System behaviors are usually specified with a standard TLA+ specification:

Init ∧ 2[N ext] vars ∧ L (II.1)
where Init is the initial predicate, Next is the next state action, vars is the tuple of all variables, L is a temporal formula that usually specifies a liveness constraint. N ext is a disjunction of multiple actions. An action is an expression containing primed and non-primed variables. The meaning of the non-primed variables is the variable's value in current state. The meaning of primed variables is the variable's value in the next state. An action expresses a transition relation.

Generally, important properties of a specification can often be expressed as invariants. Proving that a predicate P is an invariant of a specification means proving the formula

Init ∧ 2[N ext] v ⇒ 2P (II.2)
It is proved by finding an appropriate state predicate Inv and proving the following three formulae [START_REF] Lamport | Specifying systems: The TLA+ language and tools for hardware and software engineers[END_REF]:

Init ⇒ Inv (II.3) Inv ∧ [N ext] v ⇒ Inv (II.4) Inv ⇒ P (II.5)

III. SPECIFYING CONVEX HULL AND GEOMETRIC ALGORITHM

A. Specifying Convex Hull

Intuitively, the convex hull of a set of points in the plane is the shape taken by a rubber band stretched around nails pounded into the plane at each point. Formally, we give a slightly different definition from [START_REF] Pichardie | Formalizing convex hull algorithms[END_REF] in TLA+2: . Definition 3.1: Let P and S be two arrays of points, we say that S is a convex hull of P if the following predicate holds,

IsCnvxHull(S, P ) SubArray(S, P ) ∧ Len(S) ≥ 2 ∧ Dif f erent(S) ∧ Dif f erent(P ) ∧ Inside(P, S) (III.1)
where SubArray(S, P ) is an array-based alternative representation to subset, which is defined as follows SubArray(S, P ) ∀i :

∃j : i ≥ 1 ∧ i ≤ Len(S) ∧ j ≥ 1 ∧ j ≤ Len(P ) ⇒ S[i] = P [j]
Dif f erent(P ) decides whether all the points in P are different, Inside(P, S) is used to decide whether all the points of P lies inside of the polygon formed by S,

Inside(P, S) ∀i, j : i ≥ 1 ∧ i ≤ Len(S) -1 ∧ j ≥ 1 ∧ j ≤ Len(P ) ⇒ N otAtRight(S[i + 1], S[i], P [j]) ∧ N otAtRight(S[1], S[Len(S)], P [j])
where N otAtRight(A, B, C) means C does not lies to the right of the directed line from A to B. This formal definition will be used to define the postcondition of a convex hull algorithm. Obviously, the postcondition will consist of many predicates with quantifiers which is exactly what we are interested in.

B. Specifying Graham's Algorithm

The pseudo code of Graham's Scan can be found in [START_REF] O'rourke | Computational geometry[END_REF]. We specify it in TLA+2 as follows. Remark that we omit the sort operation in this specification and assume that P is already sorted, which means that all the points appear strictly counterclockwise around the point P [START_REF] Preparata | Computational geometry: an introduction[END_REF] in the order in which they appear in P . S is a stack used to store the resulting convex hull. According to the specification, the verification goal "when the algorithm terminates, the point array S is the convex hull of P " can be formalized in TLA+2 as:

P rop pc = "Done" ⇒ IsCnvxHull(S, P )
It is regarded as an invariant that should always hold for the specification. As a result, we obtain the final Verification Condition(VC) of the algorithm that corresponds to formula (II.2).

V C Init ∧ 2[N ext] vars ⇒ 2P rop (III.2)
As explained in Section II, to prove this formula, we need to find and prove an inductive loop invariant that implies this one. Finding a loop invariant is out of the scope of this paper. In the following section, we focus on how to prove a complex inductive loop invariant by the proposed decomposition rules.

IV. PROVING GEOMETRIC ALGORITHMS

In this section, we take Graham's algorithm as example to demonstrate on how to prove the correctness of a geometric algorithm. Given a geometric algorithm, we follow the paradigm described by the formulae (II.3), (II.4), (II.5) to perform the proof procedure.

A. Establishing the Loop Invariant

A standard method to obtain a candidate loop invariant is replace constants occurring in the postcondition by some variables involved in the loop and then revise the candidate loop invariant according to the failed proof attempt [START_REF] Stark | Invariant discovery via failed proof attempts[END_REF]. Inspired by the standard method, we obtain heuristically the initial candidate loop invariant IsCnvxHull(S, F irst(P, i -1)) from the postcondition IsCnvxHull(S, P ), where F irst(P, i -1) means < P [1], ..., P [i-1] >. Through multiple times of proof attempt and revision, we reach the final loop invariant.

I ∀j : j ≥ 2 ∧ j ≤ i -1 ⇒ (P [j] = S[1] ⇒ IsCnvxHull(S, F irst(P, j))) ∧ ∀j, k : j ≥ 2 ∧ j ≤ i -1 ∧ k ≥ j + 1 ∧ k ≤ i -1 ⇒ (P = S[1] ⇒ ¬Lef tOf (S[1], P [k], P [i])) ∧ i ≤ Len(P ) + 1 ∧ i = Len(P ) + 1 ⇒ S[1] = P [Len(P )] ∧ pc = "Done" ⇒ i = Len(P ) + 1 (IV.1)
According to Definition 3.1, the expansion of the formula (IV.1) consists of a number of predicates with quantifiers. This is a typical form for most geometric algorithms. In the following subsection, we introduce how to prove a loop invariant of this form.

B. Decompose Proof Obligation

According to the paradigm introduced in Section II, the proof of the formula (III.2) is transformed into the proof of the following three formulae (II.3), (II.4) and (II.5).

Init ⇒ I (IV.2) [N ext] vars ∧ I ⇒ I (IV.3) I ⇒ P rop (IV.4)
where Init and N ext are defined in the Graham's algorithm specification. The proofs of formula (IV.2) and (IV.4) are relatively intuitive, we don't bother to detail on them. We lay the emphasis on the proof of the formula (IV.3) which is much more complicated due to the existence of numbers of predicates with quantifiers. Regarding this kind of proof obligation, we propose a structure decomposition rule and a quantifier decomposition rule correspondingly for the overall proof structure and the predicates, both of which can be implemented on computer to break a monolithic proof obligation into many simpler cases. Rule 4.1: Structure Decomposition Rule(SDR). For a given a geometric algorithm, let N ext be the next state action and Inv be the candidate loop invariant. Suppose N ext has the following DN F and Inv has the following DN F and CN F respectively:

N ext A 1 ∨ A 2 ∨ ... ∨ A p Inv B 1 ∨ B 2 ∨ ... ∨ B q C 1 ∧ C 2 ∧ ... ∧ C r then the proof obligation [N ext] vars ∧ Inv ⇒ Inv is equivalent to the following formula: ∀j, k, l :j ≥ 1 ∧ j ≤ p ∧ k ≥ 1 ∧ k ≤ q ∧ l ≥ 1 ∧ l ≤ r ⇒ (A j ∧ B k ⇒ C l ) (IV.5)
According to the rule SDR, by computing the corresponding DNF and CNF of the specification and the loop invariant, a monolithic proof obligation can be broken into p × q × r simpler proof obligations, where p, q and r are decided by the depth of the expansion of the predicates. In order to control the scale of the generated subproof obligations, the expansion should be performed case by case and layer by layer according to the result of each subproof attempt.

Based on this method, for example, the DNF of Graham's algorithm specification consists of 4 conjunctions (refer to the TLA+2 specification of the algorithm in Section III-B):

1) A 1 pc = "Lbl 1 " ∧ i ≤ Len(P ) ∧ Lef tOf (S[2], S[1], P [i]) ∧ S = P ush(S, P [i]) ∧ i = i + 1 ∧ pc = "Lbl 1 " 2) A 2 pc = "Lbl 1 "∧i ≤ Len(P )∧¬Lef tOf (S[2], S[1], P [i])∧ S = P op(S) ∧ i = i ∧ pc = "Lbl 1 " 3) A 3 pc = "Lbl 1 " ∧ ¬(i ≤ Len(P )) ∧ i = i ∧ S = S ∧ pc = "Done" 4) A 4 pc = "Done" ∧ UNCHANGED vars
As for the loop invariant I, the DNF has 1 conjunctions and the CNF has 5 disjunctions in the case of non-expansion. According to the SDR rule, the number of the resulting sub-proof obligations is 20. Among of these subproof obligations, some are fairly trivial and hence can be proved automatically, some need to be further decomposed.

In the following part, we introduce another rule-quantifier decomposition rule. The idea is based on the observation that every combination of the partitions of the domains of the quantified variables reflects a different geometric situation that may require different geometric theorems to prove.

Suppose C is a formula with the following PNF:

C Γ s : s ∈ S ⇒ P ( s) (IV.6)
where

Γ (Q 1 , Q 2 , ..., Q n ), Q i ∈ {∀, ∃}, s (s 1 , s 2 , ..., s n ) Γ s Q 1 s 1 : Q 2 s 2 : ... : Q n s n S S 1 × S 2 × . . . × S n (i.e. = n i=1 S i )
Under this definition, the n-dimension space S can be represented as the union of 2 n n-dimension subspaces, which can be written as:

S = S i ∈{S i ,S * i } i=1...n ( n j=1 Sj) (IV.7)
where Si Si ∩ S i , S * i S i -Si ∩ S i .

Consider that the predicates with universal quantifiers have different decomposition paradigm from the predicates with existential quantifiers, we divide S i into two groups by the quantifier Q i , represented as:

S ∀ {i|Q i = ∀}, S ∃ {i|Q i = ∃}
Rule 4.2: Quantifier Decomposition Rule(QDR). Suppose A is an action of the algorithm specification, B is a conjunction of the DNF of the loop invariant, C is a disjunction of the CNF of the loop invariant, which has the PNF expressed as the formula (IV.6). For a given subproof obligation A ∧ B ⇒ C , it's equivalent to the following formula:

Si∈{Si,S * i } i∈S ∀ ( Sj ∈ {Sj ,S * j } j∈S ∃ A ∧ B ⇒ (Γ s : s ∈ n k=1 S k ⇒ P ( s))) (IV.8)
According to the rule QDR, a subproof obligation is decomposed into 2 |S ∀ | disjunctions.

Usually, proving every conjunct of the above formula requires different proof strategy depending on the essence of the combination of the partitions. In fact, proving the loop invariant thoroughly involves several important geometric theorems and plenty of human intervention, decomposition provides an automatic mechanism to separate the geometric situations in which different geometric theorems apply.

V. EXPERIMENT

In this proving practice of Graham's algorithm, we wrote more than 800 lines of proof codes in the integrated development environment TLA+ toolbox with the version 1.2.1. We ran it on a laptop with an Intel Core I3-370 2.4G CPU and 4G memory space. It took the computer 4.5 minutes to finish all the proofs.

VI. RELATED WORK

The application of formal method to computational geometry problem is not pervasive. In 1992, Knuth designed a counter-clockwise(CC) system by an axiomatic approach in his book Axioms and Hulls [START_REF] Knuth | Axioms and Hulls[END_REF] to formalize geometric problem. Based on the geometric axioms presented by Knuth, Pichardie and Bertot proved the correctness of an incremental algorithm and a package wrapping algorithm with the theorem prover Coq [START_REF] Pichardie | Formalizing convex hull algorithms[END_REF]. Chou [START_REF] Chou | Machine proofs in geometry[END_REF] uses algebraic method to automatically prove theorems about figures like straight lines, triangles, and circles, geometric algorithms are out of their scope.

Verifying an algorithm by verifying its loop invariant by the theory of Floyd-Hoare logic is one of the main methods. Denmat et al and Ponsini et al proposed constraint-based reasoning approach to prove the loop invariant [START_REF] Denmat | Proving or disproving likely invariants with constraint reasoning[END_REF], [START_REF] Ponsini | Automatic verification of loop invariants[END_REF]. Some heuristic methods about finding loop invariant are proposed in [START_REF] Stark | Invariant discovery via failed proof attempts[END_REF], [START_REF] Kauer | Mechanical inference of invariants for FOR-loops[END_REF], [START_REF] Ponsini | Automatic verification of loop invariants[END_REF].

VII. CONCLUSION

In this paper, we propose an approach to proving the correctness of geometric algorithms. Focusing on the loop invariants consisting of some predicates with quantifiers, we propose two proof obligation decomposition rules to transform a monolithic proof obligation into numbers of simpler subproof obligations, which can implemented to improve the automation of the proof procedure. By taking Graham's algorithm as example, we demonstrate the application of the decomposition rules. The result shows that the approach is practical in computational geometry field.

In the future work, we will try to improve the automation extent of the approach and implement it in tools.

  VARIABLES i, P, S, pc vars == << i, P, S, pc >> Init == /\i = 3/\S = <<P[2],P[1]>>/\pc = "Lbl_1" Lbl_1 == /\ pc = "Lbl_1" /\ IF (i <= Len(P)) THEN /\IF(LeftOf(S[2],S[1],P[i])) THEN /\ S' = Push(S,P[i]) /\ i' = i + 1 ELSE /\ S' = Pop(S) /\ UNCHANGED i /\ pc' = "Lbl_1" ELSE /\ pc' = "Done" /\ UNCHANGED << i, S >> /\ UNCHANGED << P >> Next == Lbl_1 \/ (pc = "Done" /\ UNCHANGED vars) Spec == Init /\ [][Next]_vars PostCondition == IsCnvxHull(S,P)

ACKNOWLEDGMENT

This research is supported in part by NSFC Programs (No.91018015,No.60811130468) and 973 Program (No.2010CB328003) of China and system@tic-Ile de France Region-CSDL.