Having access to CSD-Materials, one can predict the morphology using the Bravais, Friedel, Donnay, and Harker (*BFDH*) crystal morphology tool (Figure 5). This method implies the idea that the crystals are preferable to grow along the direction with strong intermolecular interactions. Not only simulated but also experimental morphology can be depicted and analyzed to evaluate functional groups which form a crystal face and to estimate preferable interactions which occur between this face, a solvent, or an additive. Crystal morphology prediction can be used to estimate the direction of a preferred orientation for an experimental powder XRD patterns. Photomicrographs were taken during isothermal dehydration combined with calculated BFDH morphology to help distinguish concomitant polymorphs [156], to understand the faces, and directions of solvent loss [157], while atomic force microscopy supported with BFDH calculations was used to follow phase changes *in situ* [158]. Analysis of intermolecular interactions over molecular surfaces allows the rationalization of concomitant polymorphism of substances where single crystals of two crystal forms share faces of their crystals [159]. The case of intergrowth polymorphism [160,161] referred to the existence of two solid forms within one single crystal, which is a special case of such concomitant polymorphism. On the contrary, absence of similar surfaces on dominant phases of 1:1 and 3:2 co-crystals of *p*-toulenesulfonamide and triphenylphosphine oxide allowed Croker, et al. to conclude that solvent mediated phase transition in this system occurs through dissolving of one phase and re-crystallization of the other [162]. Analysis of polymorphs, solvates, and hydrates of trospium chloride containing structures demonstrated that these compounds can be divided into three main structural groups, and the predicted and experimental crystal morphology of various forms within one structural group was found to be very similar [113]. The BFDH morphology tool was applied to rationalize silicon oil induced spontaneous phase transition in ethynyl-substituted benzamides [163] and unusual mechanical response from a crystal undergoing topochemical dimerization [164].

More complex applications of BFDH morphology tool include synthesis of desired solid forms or prediction of desired properties for previously published compounds. It was shown that self-assembled monolayers of rigid biphenyl thiols [165–167] or silioxane-based monolayers [168,169] can be used as templates to exclude concomitant polymorphysm (Figure 6a) or to obtain a metastable polymorph. First, analysis of crystal faces of the desired solid form was carried out to reveal functional groups on one of the main faces; then, the monolayer containing functional groups are able to form stable synthons with the above functional groups, which was used during crystallization (Figure 6b). Not only predesigned monolayers but also polymers [170,171] and additives (Figure 6c) [172,173] can be used as templates for crystallization. For example, β-1,4-saccharides act as templates to produce the metastable form III of paracetamol at crystallization from melt readily forming H-bonds with the (010) surface of form III [172], while the presence of sulfamides in solution promotes crystallization of γ-pyrazinamide from aqueous solutions [173]. Pyrazineamide usually nucleates from solutions in its' α form that contains molecular H-bonded dimers; amide of pyrazineamides probably forms heterodimers with sulfamide groups of templates on one of the crystal surfaces and then forms N-H...N connected chains with neighboring molecules. These chains are found only in γ-polymorph. Patel, Nguyen & Chadwick [171] not only used polymers to promote heterogeneous nucleation but also suggested to use *H-bonding propensities tool* to rank polymer surfaces towards heterogeneous nucleation of benzocaine and 1,1 -bi-2-naphthol.

**Figure 6.** Schematic representation of (**a**) concomitant polymorphism, (**b**) crystallization on self-assembled monolayers, (**c**) crystallization on additives.

Applications of some anisotropic (needle, plate) crystals are related with anisotropy of chemical bonds within such solids. Anisotropic crystals with predominated faces formed by hydrophobic groups can be used as superhydrophobic porous materials [174]. Analysis of mechanical stimulus applied along three axes of unit cells to some layered energetic materials allows us to conclude that these explosives can convert kinetic energy into layer sliding to prevent the formation of hot spots [175]. This should also be the rationale for why high-performance insensitive energetic materials can be used as desensitizers versus mechanical stimuli. Vice versa, a corrugated molecular structure of celecoxib without any slip-planes and numerous weak interactions in orthogonal directions were suggested to be the reasons for exceptionally high elasticity of its needle singe crystals [176]. Flat naphthalene diimine derivatives were found to have comparable crystal packing governed by π...π stacking, similar crystal morphology, but exhibited various mechanical flexibility attributed to close packing or interlock of their terminal alkyl chains [177]. Investigation of sulfa drug crystals [178], co-crystals of vanillin isomers [179], and amino acids [180] demonstrated that H-bonded layered structures with orthogonal distribution of strong and weak interactions attain the feasibility of cleaving a crystal along a given crystallographic plane parallel with these robust synthons.Particularly, single crystals with molecular surfaces formed by hydrophobic interactions can be applied as clean surfaces for molecular beam epitaxy. The presence of slip planes is also thought to be associated with low elastic recovery upon compression, and greater plasticity [181,182]. Crystal morphology prediction was applied to discover new gelators by focusing on scaffolds with predicted high aspect ratio crystals [183].

Thus, analysis of crystal morphology and functional groups on main crystal faces cover many areas of crystal design and are important for many areas of industry. The tools implemented to the Mercury package can be helpful in these areas, but their application seems to be limited with anisotropic crystals or solids faceted with different functional groups. In other cases, this approach can be improved by DFT calculations of the interaction effect between the growing faces and the solvent molecules, or molecular electrostatic potential on crystal faces.

#### **5. Knowlede-Based Prediction of Supramolecular Associates and Solid Forms**

#### *5.1. Polymorph Assessment*

Understanding and prediction of polymorphism—the ability of a solid material to exist in multiple crystal forms known as polymorphs—is vital for the industry as polymorphs of drugs, explosives, and pigments exhibit different properties [55]. It can be associated with various molecular conformations of flexible molecules, competing intermolecular interactions of polytopic molecules, various packing of molecules or supramolecular synthons, et cetera. While some paths to control polymorphism are described in Section 4.1, here, the approach to estimate the likelihoods of H-bonded polymorphism is described. It is of special interest for the analysis of compounds containing few donors (D) and acceptors (A) of H-bonding, since competition of donors or acceptors to take part in a hydrogen bond, and restrictions to form all highly likely bonds simultaneously may cause the appearance of polymorphs [184]. Taking into account that various polymorphs exhibit different properties including solubility, bioavailiability, tabletability, et cetera, the pharmaceutical industry has prompted the development of computational tools that are able to predict the likelihoods of H-bonded polymorphs. Thus, a methodology has been developed [185] to estimate the likelihood of H-bond formation for each pair of donors and acceptors in a molecule taking the environment of the D and A groups into account.

The *H-bond propensities tool* uses only a 2D molecular formula (Figure 7) and involves four stages: Data sampling, model fitting, model validation, and target assessment. For each functional group, it is either automatically assigned from the 2D formula, taken from the functional group library, or sketched manually by using a fitting data set, which is generated by loading in an existing set of structural data, or from the CSD. The H-bonds within the fitting data set are identified to collate statistics and descriptors of H-bonding formation. The logistic analysis is performed to generate a statistical model to determine the likelihood of H-bond formation; the model should be analyzed to include for the final model only significant variables and be good enough to proceed (area under ROC curve (receiver operating curve) above 0.8 indicates good discrimination). The results are ranked by propensity to allow inspection of the most and the least likely D-A pairs. Based on this knowledge-based approach, a structure that adopts the highest propensity H-bonds displays a low likelihood of appearance of a more stable polymorph, while for solid forms with lower propensities of H-bonds, a high risk of polymorphism caused by H-bonding is expected. For example, the most likely interaction for a "molecule" depicted in Figure 7 includes a donor atom of Group 3 and an acceptor atom of Group 1. However, in solid {Refcode} the second likely hydrogen bond between groups 1 and 2 is observed, thus, indicating possibility of appearance of another polymorph with the most likely hydrogen bond.

**Figure 7.** Flowchart for the methodology to predict the likelihood of hydrogen bond formation.

The tool has proved to be highly valuable in estimating the relative stability of known crystal forms [186]. Particularly, the kinetically more favored **Form I** of ritonavir displays statistically unlikely hydroxy-thiazoyl and ureido-ureido interactions but exhibits more favorable conformation of the carbamate moiety, while **Form II** realizes stronger hydrogen bonds as estimated in Ref. [186]. Two X-rayed polymorphs of 4-aminobenzoic acid realize different H-bonding, and the more stable form realizes more likely sets of H-bonding [187]. The approach was also successfully used to rationalize polymorphism for N2-(indol-3-acetyl)-L-aspargin [186], lamotrigine [188], and crizotinib [189,190] (see Figure 8 for chemical formulas of these compounds). Not only strong donors and acceptors of

H-bonding can be used for analysis; a competition of hydroxide O-H and ethynyl C−C≡C−H groups to form hydrogen bonds [191] and the role of CHCl3 and CH2Cl2 molecules in crystal packing [192] were demonstrated using this tool.

**Figure 8.** Schematic representation of compounds analyzed using H-bonding propensities tool.

Known forms of paracetamol both exhibit H-bonding between OH groups, that is a less likely donor to form a bond than an amide fragment. It was assumed that there was a conflict between a good donor group to form a bond disfavoring a poor acceptor to form bonding, which could be overcome by co-crystal formation with molecules containing the same D and A groups [187]. Indeed, co-crystals of paracetamol with diamino- or bis(4-pyridyl)-containing co-formers contain molecules of paracetamol connected by more likely amide/hydroxo or amide/amide pairs of interactions [193]. A good correlation between H-bonding charge analysis and H-bond propensities was also demonstrated for crizotinib [189], and heterocycle-1-carbohydrazoneamides [194].

Note, that the high probability for exhibiting polymorphism does not necessarily allow one to obtain novel polymorphs and that different polymorphs can realize the same H-bonds [195]. For example, probenecid was predicted to realize various H-bonded polymorphs, and differential scanning calorimetry indeed revealed three polymorphs, but XRD experiments showed that all of them realized the same hydrogen bonding [196]. For bufexamac, the propensity tool suggests the possible existence of three crystal forms, while Nauha and Bernstein found two polymorphs with similar H-bonds [197]. Two experimentally observed forms of meglumine realize slightly different H-bonding, which, however, are absent on putative structure landscape as these contain highly unlikely bifurcate acceptor OH groups [197]. For axitinib [189], five anhydrous polymorphs realize the most likely H-bonding, thus, the approach can not distinguish the most stable among them.

To sum up, the propensity tool neither predicts polymorphism, nor guarantees that all sets of theoretically possible H-bonding combinations can be obtained or guidelines can be given on how to obtain any of these polymorphs. It does not describe inclinations of a molecule to form concomitant polymorphs, conformational polymorphs, or packing polymorphs, or solid forms with Z' >1. Instead, it indicates the possibility of solid forms to organize various H-bonded architectures and provides some guidelines on the amount of affords that one can spend in the experimental search for various crystal forms. Abramov also notes such a limitation of the approach as incapability to distinguish more stable

polymorph among crystal forms with similar H-bonded networks, which however can be overcome with additional charge-density analysis of known polymorphs [189]. At the same time, its application is not limited by monomolecular systems. Note, that any desired functional groups, additional co-formers, and solvent molecules can be included in the statistical model. Thus, although this method is not able to predict the ratio of components in a co-crystal, it can be applied for co-crystal design.

#### *5.2. Co-Crystals Design*

Co-crystals are solids that consist of two and more components (co-formers) that form a unique crystalline structure having unique properties. Thus, co-crystallization is an approach to optimize the physical properties of solid materials. For example, pharmaceutical co-crystals can be obtained to modulate dissolution rate and physical stability of drugs [198,199]. Schultheiss and Neuman reported [200] the effect of coformers on melting points, stability towards humidity, solubility, bioavailiability, and some other properties of pharmaceutical co-crystals. Karki, et al. [201] demonstrated that the tabletability of pharmaceuticals can be tuned up on the example of paracetamol co-crystals with various co-formers. The dependence of mechanical properties of acid...amid based co-crystals was studied by Saha and Desiraju [202]. The potential of co-crystallization for tuning the properties of the energetic materials [175,203,204], optical materials [205,206], and for the food industry [207] have also been demonstrated in recent papers. In contrast with synthetic routes focused on the synthesis of covalent derivatives containing functional groups that affect the desired property, co-crystallization became of interest as a 'greener' process, frequently free of toxic ingredients and by-products. Pharmaceutical compounds in this context become of particular interest as these molecules contain functional groups that can be involved in molecular recognition of biomolecules [208]. The main task of crystal chemistry in this field is to estimate which co-formers, if any, will form a co-crystal with a given molecule prior to screening and investigation of properties of the solid obtained.

One group of methods available are quantum chemical calculations from simple energy minimization to full structure prediction [209]. Another group of methods is based on data knowledge about the likely homo and heterosynthons. Since the works of Etter [210], it is known that (i) all good donors and acceptors tend to take part in hydrogen bonding, (ii) the strongest donors tend to interact with the strongest acceptors. In this context, it becomes essential for a researcher to estimate which donors and acceptors are "good" and "the strongest". The *H-bonding propensities* tool allows ranking of various donors and acceptors without any quantum-chemical calculations (Figure 9).

**Figure 9.** Flowchart for co-crystal design strategies based on synthon competition.

*H-bonding propensities* tool for uncharged and charged molecules of pyrimethamine and dicarboxylic acids, as well as halogenated aromatic compounds (Figure 10a), was successful in predicting formation and non-formation of adducts [211]. A similar approach was later applied to pyrimethamine with some other drug molecules also taking solvent molecules into account [212]. Note, that proton transfer accompanied with transformation of a co-crystal to a salt was rationalized in terms of their Δp*K*a values [213,214], where ΔpKa = pKa[protonated base] − pKa[acid] >2 or 3 was found for salts, negative Δp*K*a was characteristic of co-crystals, while an intermediate situation may result in both depending on stable supramolecular synthons [215]. Hydrogen-bond propensities gave similar trends as calculated bond energies for possible synthons in thiazole amides (Figure 10b), and co-crystallization of six thiazole amides with 20 different carboxylic acids demonstrated effectiveness of complimentary energy-based and data–knowledge predictions for the prediction of likely synthons [216]. The knowledge-based approach to co-crystal design was successfully applied to select pyrazinoids and pyridinoids as prospective co-formers for a diuretic drug hydrochlorothiazide [217] to tune its solubility. The amide-pseudoamide synthon was shown to be more stable than the formation of two dimers of the theophylline molecule, that allowed to explain the formation of theophylline: Amide co-crystals [218]. At the same time, theophylline co-crystals with fluorobenzoic acids demonstrate that it is still difficult to predict co-crystallization of small rigid molecules using this method [219]. Similarly, preferable co-crystallization of only one of three tautomers of 2-amino-6-methyl-1,4-dihydropyrimidin-4-one with carboxylic acids was explained using this approach [220]. Potential of co-crystal design for molecules with limited H-bonding functionalities was demonstrated on the example of propyphenazone [221]: NH2, OH, and CO2H functionalities were found to be the most likely groups to interact with O=C group, and eight novel co-crystals were obtained using co-formers containing OH or/and CO2H groups. Dicarboxylic acids [222] or anions [223,224] can be applied to co-crystallize two molecules, which do not form co-crystals.

**Figure 10.** Examples of homo/heterosynthons analyzed for (**a**) pyrimethamine [211,212], and (**b**) thiazole-amides [216].

Note, that the analysis of possible heterosynthons is not limited to co-crystals. It can be utilized to investigate inclusion compounds [225] and mixtures of polymers and solutes [171]. Besides, the ranking of heterosynthons on their relative strength can not only help in estimating the most likely binary co-crystals. This method has great potential in the synthesis of multicomponent co-crystals based on polytopic co-formers. Thus, H-bonded heterosynthons are able to interact with the third co-former through halogen bonds [226–228]. π-Stacking in conjunction with hydrogen bonding was used for synthesis of ternary co-crystals [229,230]. Combination of H-bonded synthons and stacking interactions between donor and acceptor planar molecules allowed Desiraju and co-workers to synthesize quaternary and even quintinary co-crystals [231,232]. Partial substitution of some co-formers in quaternary systems with their shape-size analogs even allowed for obtaining six component solids [233].

Although some examples of co-crystals design based on H-bonding propensities tool are given above, successful prediction of inclination of some molecules to form co-crystals obviously needs analysis of some other molecular descriptors. First, the possibility of a solvent to take part in H-bonding should be taken into account [211,212,234]. A systematic study of co-crystallization of paracetamol with H-bond acceptors [235] and donors [236] demonstrated that paracetamol molecules in these crystals are linked via either OH··· O=C or NH··· O=C interactions, depending on the presence or absence of substituent groups on the molecule of the second co-crystal former. A machine learning algorithm trained out on a set of paracetamol co-crystal experiments using more than 190 molecular descriptors for each co-former allowed to predict 9 of the 13 experimentally obtained co-crystals within the top 11 suggestions [237]. Unfortunately, this method requires a large amount of experimental work to train the model, and the model obtained should be applied only to molecules with similar molecular descriptors as those for an attested molecule. Fabian has proposed a methodology that involves using size and shape complementarity as the primary driver for co-crystal formation [238]. Analysis of 131 molecular descriptors for 710 co-crystal partners suggested that co-crystals were more likely to form between molecules of a similar size and shape and with a similar polarity of co-crystal formers, thus both these descriptors were recently included to the *Co-Crystal Design tool*.

Effectiveness of such modification can be demonstrated on the example of 1,2,4-thiadiazole derivative co-crystallization with gallic and vanillic acids [239]. While *H-bonding propensities tool* gives almost equal probability of occurrence of homo- and heterosynthons, molecular complementarity tool indicates that co-crystals in these systems should form; and these co-crystals were experimentally obtained. Karki, et al. [240] demonstrated that synthon analysis combined with Fabian's methodology was effective in prediction of possible co-formers for artemisinin. A series of co-crystals of sulfamethoxazole [241] and leflunomide [242] were synthesized using this approach. Prediction of possible H-bonded motifs between tyraminium cations and violurate anions also successfully predicted many of the bimolecular synthons experimentally observed in tyraminium violurate polymorphs and hydrates, but also demonstrated that none of the trimolecular synthons were predicted [243]. Note, that tri- and tetra-molecular synthons are not something unusual in co-crystals (see, for example, tetramolecular associates found in thiazole-amides, Figure 10), thus any predictions should also include such polymolecular associates into account.

To sum up, *H-bonding Propensities* and *Co-Crystal Design* tools allow one to simplify co-crystal screening. CSD analysis helps to select complementary functional groups to form heterosynthons more likely to form than homosynthons. Presence of false positives in experimental screening then may be accounted for other factors. For example, although hydrogen bonds are thought to be stronger than halogen bonding, one can undergo a competitive co-crystallization between H-bonded and halogen bonded synthons based on the polarity of a solvent used [244]. It was shown that in polar solvents 1,2-bis(pyrid-4-yl)ethane forms co-crystals via halogen bonds, and in nonpolar - via H-bonds. Of more concern with regards to co-former screening is the possibility of false negatives. If the results of this type of analysis are used to narrow down the number of screening experiments performed, any possible co-crystals that are incorrectly marked as unlikely to form would be missed. Besides quantum calculations of molecular dimers or logistic models trained on large datasets of co-crystals, development of shape-size molecular complementarity tools for the prediction of possible co-formers seem to be very promising. Note, that the CSD-Crossminer package developed for sophisticated search of analogues of drugs is now available (Section 3.4). It utilizes the idea that molecules with similar disposition of functional groups and hydrophobic fragments may act similarly with proteins and replace each other. The idea of substitution of one of the co-crystal formers with its size-shape analog can be realized using this program.

#### *5.3. Full Interaction Maps*

The *Full Interaction Maps* (FIM) tool [245] implemented within the Mercury package visualizes the likelihood of a synthon appearance and corresponding geometry variation between functional groups of a molecule under consideration, and a probe functional group. For each functional group of a molecule, 3D scatterplots of CSD contact searches with a chosen probe, and the functional groups are generated and converted into scaled density maps. These density maps are then combined for the whole molecule, taking the environmental effects of combinative factors and steric exclusion to account. The 2019 version of the CSD-Enterprise contains (i) RNH3, uncharged and charged NH nitrogen atoms as probe functional groups of donors of hydrogen bonds, (ii) various oxygen atoms (of a carbonyl and alcohol group or a water molecules) as a probe of hydrogen bond acceptors, (iii) methyl and aromatic carbon atoms as a probe of stacking interactions. The CSD-Enterprise version of 2018 allowed the use of C-I and C-Br interactions as a probe of (iv) halogen interactions; and the next year these interactions were extended with C-F and C-Cl probes, to evaluate differences between various halogen atoms. For each functional group, the FIM distribution is similar with *IsoStar* scatterplots of corresponding

functional groups, but the FIM for a whole molecule reflects the fact that some of functional groups are better donors and acceptors of, for example, H-bonding, than the others. On these maps, this difference is expressed in the color (and the coordinates of the most likely positions of the donor and acceptor sites can be additionally calculated and depicted), while corresponding quantitative values can be estimated using the *H-bonding Propensities* tool (Section 5.1). Besides, these maps are very sensitive to the steric exclusion from other molecular species and would differ for a molecule in optimized and experimental geometries, and those obtained with CSD-Conformer Generator.

The FIMs probed with H-bond donors (blue) and acceptors (red) and supramolecular synthons of the *Pbcn* polymorph of chalcone {BZYACO01}and 5-(4-bromophenyl)-1-phenylpent-2-en-4-yn-1-one {YILNAK} are depicted in Figure 11. For both molecules' blue regions near the ketone oxygen atom reflect its inclination to form bifurcate bonding due to the presence of two lone pairs on the oxygen atoms. The most expected positions of acceptors of hydrogen bonding for the chalcone are situated on the opposite sites of the molecule, and correlate well with observed supramolecular synthons in its polymorphs—a dimer depicted in Figure 11, a {BZYACO01}and a head-to-tail chain {BZYACO04} similar with that of 5-(4-bromophenyl)-1-phenylpent-2-en-4-yn-1-one. Elongation of the π-conjugated chain with a triple bond makes the appearance of such chains more abundant than dimer occurrence (Figure 11b), that was experimentally confirmed for a series of 1,5-diarylpent-2-en-4-yn-1-ones [246,247].

**Figure 11.** Interaction maps and supramolecular synthons in the crystal structures of (**a**) chalcone {BZYACO01}, (**b**) (E)- 5-(4-bromophenyl)-1-phenylpent-2-en-4-yn-1-one {YILNAK}.

This tool has been applied mainly for the analysis of strong hydrogen bonds [54,170,190,194,248–251] in accord with the first test functional groups suggested within this tool, and the first paper published to describe it. Mutual disposition of hot spots can be used to find functional groups of the same molecule [245], co-formers [243,249,250], active sites in a binding pocket of a macromolecule [252], or surface inhibitors [170], which match a given pattern of interaction preferences. Analysis of polymorphs of some drugs demonstrated that disposition of H-bond donors and acceptors close to the hotspots are indicative for the more stable polymorph even if a molecule realizes less likely conformation [190]. However, in the case of polymorphism induced by weak intermolecular interaction reorganization [253] or high-pressure [254], it has lower predictive ability. This can easily be understood if we keep in mind that data for the FIMs are plotted based on structural information derived at atmospheric pressure.

At the same time, the test carbonyl oxygen and uncharged NH groups can be successfully applied to reveal positions of acceptors and donor of C-H...O [170,255–257] and C-H...N [258] bonding. It was demonstrated in Refs. [255,257] that the FIMs can be applied to estimate likely C-H...O bonded motifs for chalcones, polyenones, pentenynones and cyclic ketones with vinylacetylene fragments. At the same time, it was demonstrated that C=O and C-Br groups compete with each other for the most acidic hydrogen atoms, thus that FIM predictions of synthons based on C-H...Br bonding became less reliable at the presence of C-H...O=C bonding [257]. At the absence of carbonyl groups the FIMs for the N-salicilidenanylines probed with C-I and C-F groups successfully predict all of the C-I...N and many of the C-F...H-C interactions in co-crystals of N-salicilidenanylines with perhalogenated co-formers [259]. Mugheirbi and Tajber analyzed the FIM hotspots around the itraconazole molecule to understand the molecular environment in the mesophase [260]. The itraconazole lacks any donors of strong H-donors but has a number of competing acceptors of H-bonding; thus, analysis of the FIMs and FTIR spectra allowed them to reveal the most ordered of mesophases, and to propose that the greatest mobility of this molecule is associated with movement of the triazoline ring. Analysis of FIMs in two polymorphs of the dinuclear Co(II)-Shiff base complex revealed the most acidic H(C) atoms (similar in both polymorphs), which take part either in C-H...O bonding in the triclinic polymorph, or in the C-H...π bonding in the monoclinic polymorph, or even to not take part in any prominent intermolecular bonding [256].

Note, that unsatisfied strong acceptors observed in a crystal structure solved from powder XRD data may be indicative of missed water/solvent [261,262]. In this case coordinates of the most expected position of a water molecule can be used in the refinement. Worth noting, that *IsoStar* and FIMs deal only with intermolecular interactions but a similar idea could be used for investigation of metal-ligand bonding in polytopic ligands. This could be helpful for understanding the factors that govern linkage isomerism, and to estimate the most likely coordination mode of the most widespread ligands.

#### *5.4. CSD-Conformer Generator*

Representation of a molecule in the three-dimensional space finds numerous applications in structure solution from powder diffraction data in real space, crystal structure prediction including the formation of co-crystals, protein–ligand docking, and others. Thus, the *CSD-Conformer Generator* tool included in the CSD-Enterprise as an approach to a fast generation of plausible molecular conformations by using geometric distributions derived from the CSD. First, for the input 3D molecular model with all hydrogen atoms present, all bond lengths and angles are minimized based on corresponding average values. Then, the molecule is partitioned into rotamers, and rotamer libraries and ring template libraries are used to generate a conformer tree with unforbidden and preferred rotatable bond geometries and ring geometries. A final set of conformers is clustered according to conformer similarity. Each conformer is locally optimized in torsion space. It was demonstrated that this tool reproduced well conformations of a number of molecules observed in the CSD and the Protein Data Bank [263,264]. Some discrepancies between the predicted and experimental geometries occurred for the ligand: Macromolecule complexes in unusual conformational space, rare rotamer examples, uncertain bond types, and some other cases, nevertheless, theoretical configuration space represents well experimental data [265,266] or configurations obtained by the Molecular Operating Environment's Low Mode Molecular Dynamics module [267]. It was demonstrated that the CSD-Conformer Generator combined with ab initio [268] or DFTB3-D3 [269] calculations can be used for crystal structure prediction of some flexible pharmaceuticals. Such combination reproduced well the molecular geometry and crystal parameters, although did not provide sufficiently accurate energy ranking.

#### **6. Hydrate**/**Solvate Analyser**

As it was described above, a solvent can affect crystal morphology and sometimes determines the polymorph, but it also can be built into the crystal structure to form a hydrate or solvate. Since water is nature's solvent, non-toxic and widespread, the scientific community is interested in understanding water assembles in liquids and solids. In particular, medium-sized and large water clusters are important for biology, since they act as surrounding and solvating solutes for biologically active molecules, fill discrete voids and channels in molecular and supramolecular assemblies including their reactive sites and interpenetrate into the interfacial region of hydrophobic surfaces [270–273]. Besides, water associates are involved in dynamic processes such as proton transport, material or molecular folding, de/resolvation, and others [274,275], and even affect mechanical properties of single crystals [276]. The behavior of liquid and solid water including ice, clathrates, and ice-like systems are

determined by the disorder of the hydrogen atoms, and the binding energy of different configurations as well as some other properties can also differ [277–279]. At last, but not least different solvents may stabilize different forms of a molecule (the neutral or zwitterionic forms) [280,281]. In other words, both (i) connectivity of a solvent molecule in a crystal, inclination to form a given associate, and (ii) dimensionality and unit cell volume that goes to solvents are of interest for biology, crystallography, and material chemistry. All corresponding algorithms are realized within the *Hydrate* and *Solvate Analyser* tools of the CSD-Material module. The *Hydrate Analyser* provides information about the water H-bond geometry and motifs (one of any of 10 most common motifs described by Ref. [282]) detected, information on the volume occupied with water molecules, and display the water space and water interaction maps. The *Solvate Analyser* provides similar information about simple and mixed solvates, co-formers and ions.

#### *6.1. Analysis of the Local Connectivity of a Solvent*

Analysis of solvate connectivity and the most abundant motifs (if any) formed by solvent molecules is similar to that described above for the investigation of other synthons and supramolecular associates. The practical meaning of results obtained for material chemistry, biochemistry, and the pharmaceutical industry are based on the fact that similar hydrate architectures were found in crystals of both inorganic, organic, coordination, and even macromolecular compounds [40,41]. H-bonding between a complex and an outer sphere solvent can lead to additional quenching in luminescent materials that are unfavorable for the task-specific design of optical materials [283,284]. Comparison of clathrate, turbolato-clathrate, and non-clathrate hydrates revealed that water associates in the latter can be regarded as fragments of clathrate hydrates [285]. Analysis of hydrate architectures in crystal structures of known kosmotropic and chaotropic agents revealed that the kosmotropes tend to take part in H-bonding with hydrates, while chaotropes in crystals tend to from clathrate hydrate-like structures [286]. Kosmotropic agents are those able, by the ordering of the water structure in solution, to enhance intermolecular interactions within protein molecules, thus preventing denaturation, while chaotropes act oppositely. Hence, the property of a molecule to crystallize with water giving full enclathration was directly associated with the ability to act as chaotrope.

Most frequently water serves as a donor of two hydrogen bonds and an acceptor of one hydrogen bond, although three other motifs are also relatively common (Figure 12a [282,287,288]). The dual nature of this molecule in respect to H-bonding allowed proposing imbalance in the number of donor and acceptor groups of a polytopic molecule to be the reason for form hydrate appearance [289]. Although, it is in accordance with Etter's rule, which states that "all good proton donors and acceptors are used in hydrogen bonding" [210], Infantes, Fábián, and Motherwelldemonstrated that the imbalance of donor and acceptor groups, in fact, does not affect inclination of a molecule to form hydrates [288]. Instead, it more readily interacts with unsatisfied acceptors of H-bonding, especially R2PO2 - , Cl- , C-NH3 + groups, acts as a bridge between unsatisfied donors and acceptors, and occupies the free volume, especially for chiral molecules [290]. The strong imbalance between the number of donors and acceptor groups between tetrasulfonate-functionalized rigid anions and planar polyamino-containing cations indeed afforded their crystallization as hydrates, where water molecules act both as bridges between cations and anions, and clusters incorporated within cavities and channels of H-bonded networks [251]. Nevertheless, typically it is still hardly possible to predict if a compound will form a hydrate or not. Zaworotko and coworkers revealed that not only molecules with unsatisfied donors but also acceptors of H-bonding also readily form hydrates, but the reason of CSD statistics, in this case, might underestimate inclination of a molecule to form a hydrate, probably, because the most effective pathways to hydrate formation were slurrying in water and exposure to humidity [291]. The role of MeOH [292], DMSO [293,294], and CHCl3 and CH2Cl2[192] solvates in crystals was also investigated.

**Figure 12.** H-bonding patterns in (**a**) organic hydrates [288] and (**b**) organic peroxosolvates [116] (X denotes any donor or acceptor atom).

Comparison of motifs formed by water molecules in organic hydrates and by H2O2 in organic peroxosolvates gave a clue to the understanding of their isomorphysm [116,295]. It was found that H2O2 always forms two H-bonds as a proton donor, and up to four bonds as an acceptor of H-bonding (Figure 12b). Thus, only three of the peroxosolvate H-bonding patterns are similar to those observed in crystals of hydrates. As isomorphous substitution may occur only in crystals with similar H-bonding motifs, peroxosolvates with DD, DDA, and DDAA motifs connectivity require high H2O2 concentrations for synthesis, are sensitive to humidity, and as a result are of limited practical meaning. Instead, peroxosolvates involved in at least five hydrogen bonds and strong H-bonding with acceptors are not inclined to take part in isomorphous substitution with hydrates.

Not only water connectivity, but also the main motifs formed by water molecules, or supramolecular architectures including water molecules as well as suggested nomenclature were described [67]. Combined together, these data have practical meaning for prediction of coordinates of missing water molecules for solving and refinement of crystal structures from powder X-Ray diffraction data [261], as well as for molecular docking with GOLD (Genetic Optimization for Ligand Docking), part of the CSD-Discovery module.

#### *6.2. Analysis of Solvent Associates*

Thus, besides analysis and classification of H-bonded motifs, visualization of associates, and calculation of the unit cell volume occupied with water/solvent molecules (both coordinated and not) is available via *Hydrate*/*Solvate Analyser* tool. On Figure 13, a DMSO and water space in the structure of bosutinib DMSO solvate trihydrate {ABECES} are shown in red and blue, respectively. It is clearly seen that both solvents form isolated clusters, while in the structure of tetrabutylammonium fluoride clathrate hydrate {CIPRAV} H-bonded water molecules form a 3D architecture with tetrabutylammomium cations situated in pores of this framework (Figure 13b). Visualization is similar to void analysis, but there is no need to remove any molecules to carry out calculations. Moreover, the *Solvate analyser* tool can be used to identify not only solvents, but also co-formers and ions, disordered fragments, or even hydrophobic/hydrophilic fragmentset cetera, as well as to visualize their packing and check H-bonding involving the identified components.

Analysis of the voids within previously published copper-containing coordination polymers allowed Inokuma et al. to find a crystalline sponge with desired pore size and affinity to guest molecules and demonstrated that a huge number of previously published coordination polymers can be used to capture guest molecules [296]. Anisotropic motifs formed by solvent molecules affect mechanical properties of compounds, and dehydration mechanism. For example, presence of infinite H-bonded water chains is typically associated with the reversible dehydration process and similarity between crystal structures of the hydrated and anhydrous forms of a solid, like caffeine [297], carbamazepin [104], aspartame [298], and Shiff-base Ni(II) complex [299]. Such channels can be not only dehydrated but also substituted with other solvents [300,301]. The solid dehydration with

prominent atomic movement traced with single-crystal or powder X-Ray diffraction can shed light on a more complex mechanism of water loss [302–306]. Various mobility of the water clusters incorporated within the pores and H-bonding waters that structure these pores was demonstrated in Refs. [251,307]. Liu et al. demonstrated that anhydrous uric acid and its' dihydrate exhibited various mechanical properties of their single crystals despite the similarity of layered motifs formed by the acid [276]. 2D water associates situated between acid layers perpendicular with single crystal main faces break at indentation with dehydration at crystal surface and make this solid softer than the anhydrous form.

**Figure 13.** The structure of (**a**) bosutinib dimethyl sulfoxide (DMSO) solvate trihydrate {ABECES} and (**b**) tetrabutylammonium fluoride clathrate hydrate {CIPRAV} with solvate analyser display of DMSO space shown in red and water space in light blue.

#### **7. Additional Software Compatible with the CSD-Materials Analysis of Structure– Property Networks**

#### *7.1. CSD Python API*

The *CSD Python API* is now part of the CSD-Enterprise automatically installed as part of the CSD installation. This allows one to run supplied or user-written Python scripts for the loaded structure or a set of structures. CCDC Python-Built-In scripts include at the moment:

Some other examples of Python scripts applied for investigation of crystal structures include removing solvent molecules from a MOF subset extracted from the CSD to perform analysis of geometry and properties of nondisordered porous compounds [309]; authors realized the possibility to remove by request both uncoordinated and coordinated solvent molecules from coordination polymers. Dolinar et al. wrote a Python script to perform the CSD search of chiral monomolecular compounds [310]. Bryant, Maloney, and Sykes investigated a number of polymorphs and co-crystals with previously reported tabletabilities to analyze how different factors govern tabletability of organic crystals, and to suggest a programmatic method to calculate corresponding crystal descriptors [311]. Their investigation of the most likely slip planes, their interpenetration on connection by weaker interactions, presence of other slip planes, and automatic assessment of H-bond dimensionalities revealed that the degree of separation, followed by the presence of H-bonds between the layers, and finally d-spacing between potential slip planesare the most important descriptors affecting the tabletability. The corresponding CSD Python API protocol that can be applied for the prediction of mechanical properties of organic crystals was also published [311]. Miklitz and Jelfs reported a Python script for structural analysis of porous materials (organic and coordination cages) that estimates the cavity diameter, number of windows, and their diameter [312]. It uses Cartesian coordinates and atom types for input, thus, can be in principle applied also to selected fragments of framework materials. Besides codes to extract and analyze data stored within the CSD, for text mining of selected properties within a large set of papers was also published [313]. The algorithm uses manuscript html files as input for extraction of surface areas and pore volumes of metal–organic frameworks with at least 73% accuracy and can be used for investigation of structure–property networks.

Thus, the Python API gives a convenient path to combine and modify the CSD requests with CSD-Materials, CSD-Discovery tools or external software for ones' purposes. The next Section gives some examples of the combination of the CSD-Materials tools with external crystallographic software.

#### *7.2. Combination of Various Tools and Algorythms*

Some combinations of the CSD-Materials tools with external software were mentioned in previous Sections, for example, the combination of the *CSD Conformer Generator* with periodic calculations to predict crystal structures, or that of the *CSD Co-Crystal Former* and the *H-bond propensities* with DFT calculations to better understand competing for intermolecular interactions. Let us mention also some other examples, where the *CSD-Materials* tools were used in unusual way or as the tool to solve only part of a scientific problem.

Crystal packing of organic single-component compounds with Z' = 1 were compared to evaluate the effect of chemical transformation on crystal packing and isostructurality using a combination of the CSD Python API scripts, external software, and the *CSD Crystal Packing Similarity* tool [85]. For this purpose, 15,5543 organic compounds were extracted from the CSD based on their composition and data quality. Then, a freely available algorythm published in 2010 [314] was used to cleave up to three acyclic single bonds between functional groups as demonstrated in Figure 14. Comparison of the *Crystal Packing Similarity* between molecular clusters formed by the "Value" molecular parts, allows us to reveal the most frequent terminal transformations of the "Key" groups and to evaluate their effect of the packing and isostructurality. For example, substitution of the methyl group to the chloro-, azide- or iodo- groups gives isostructural pairs in nearly 30% of cases.

**Figure 14.** (**a**) Flowchart of the algorythm of bond cleavage reported at Refs. [85,314] and (**b**) some fragments formed from single, double, and triple cuts of lamotrigine.

Recently possibility of prediction of H-bonded motifs using a topological approach was demonstrated for some families of organic molecules with a large number of X-rayed representatives [315,316]. This method schematically represented in Figure 15 utilizes the idea that (i) only a limited number of architectures can be formed from a particular building block with a given connectivity, (ii) overall topology of a network depends on connectivity of building blocks, and (iii) propensities of various networks appearance are not equal. While relations between connectivity of an organic molecule and resulting network and corresponding propensities of their appearance have already been published [315,317], this approach does not allow evaluating the most abundant molecular connectivity of insufficiently studied families of organic compounds. Recently, maps of electrostatic potential, the *Full Interaction Maps*, and the *H-bonding Propensities* were compared as the tools to overcome this problem for heterocycle-1-carbohydrazoneamides [194]. Knowledge-based approaches demonstrated competition between various nitrogen atoms to act as acceptors of H-bonding with the donor -NH2 group and allowed evaluating local molecular connectivity. Experimental crystal structures realize one of highly abundant theoretically predicted topologies of an H-bonded network.

**Figure 15.** Flowchart of the algorithm to estimate H-bonding architectures of organic molecules.

Tilbury, et al. suggested to combine the *H-bond propensity* based prediction with COSMO-RS theory [318,319] to predict drug substance hydrate formation [320]. Hydrate formation probability is estimated based on propensities calculated for each donor/acceptor groups in drug...drug or drug...water pairs. The most favorable hydrate formation has the maximal difference between these values. Although within this model relative strength of different functional groups to interact

with water is displayed, this model does not take steric effects and intramolecular interactions into account. Additional quantum mechanics approach to eliminate donor and acceptor groups that are unable to form intermolecular interactions and to calculate the most favorable molecular and hydrate conformations allowed to improve predictive fidelity of calculations.

Chandra, et al. [225] applied CSD to enhance the solubility of telmisartan (TEL), a low soluble antihypertensive drug. First, the possibility of TEL to form likely intermolecular interactions with sulfobutylether beta-cyclodextrin was confirmed using H-bonding propensities tool. Then, the best docking pose of an inclusion complex of TEL in cyclodextrine was found with GOLD. The predicted binding mode was in line with the experimental spectra for the inclusion complex, which demonstrated enhanced solubility and dissolution rate of TEL.

Correlation between solubility and various parameters describing dimethyl sulfoxide (DMSO) role in crystals structures was carried out by Spitery, et al. [294]. They analyzed various parameters of hydrogen bonding of DMSO molecules as obtained from crystal structures reported in the CSD with solubility estimated using Chem3D [321], and found a negative correlation between the number of interactions the solvent is involved in and solubility of a solid.

#### **8. Conclusions**

The Cambridge Crystallographic Data Center holds a unique position of a curator of the Cambridge Crystallographic Database that now contains more than 10,000,00 crystal structures and a software developer for the analysis of collected data. This software allows the visualization of complex statistical data in a click-of-a-button manner (*Full Interaction Maps, Conformer Generator, Mogul analysis, Solvate Analyzer, BFDH predictions*), and investigation of local connectivity of functional groups and small molecules (*H-Bond Propensities, Co-Crystal Former, Motif Search*) that became possible due to the network between chemical diagrams and automatically calculating interatomic and intermolecular distances. Despite the complexity of the utilized algorithms, big data analysis, and their sophisticated processing, the software remains visible, user-friendly, and available for each CSD user. Thus, its' visibility and flexibility combined with exhaustive free documentation (see user guides, how-to-video, tutorials, and other materials at [322]) provide the basis for successful application in crystal engineering, material chemistry, and chemoinformatics.The diversity of CSD-supported studies described in this review indicate the great potential of the knowledge-based software for future development in these and associated fields, thus, the development of links with other molecular and crystallographic databases and repositories of properties would be advantageous.

The combination of the CSD tools with external software similarly with UNI and MOPAC calculations incorporated within the CSD-Enterprise seems to be very prospective for future development. DFT calculations, especially periodic ones, could provide further insight into the field of coordination polymers with numerous applications in gas storage and separation, catalysis, and spintronics. Implementation of graph theory to polytopic molecules and ligands could expand data about the local (molecular) connectivity to knowledge about possible and most abundant coordination and H-bonded architectures. The future of this field still possesses many challenges, but there is no denying that the value of the Cambridge Structural Database and associated software for the knowledge-based analysis in the filed of composition–structure–properties studies can not be overestimated.

**Funding:** This work was supported by the Ministry of Science and Higher Education of the Russian Federation. **Acknowledgments:** The author is grateful to A. A. Korlyukov for fruitful discussion of the text.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**

© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Crystals* Editorial Office E-mail<EMAIL_ADDRESS>www.mdpi.com/journal/crystals

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18


*2.6. ompF and ompC Expression within Y. pseudotuberculosis Population*

Current experimental evidence suggests that bacterial cultures are constituted of heterogeneous subpopulations [55] with differential responses to environmental stresses [56,57]. For instance, Sánchez-Romero and Casadesús have observed the population heterogeneity of *S. enterica* cells in the expression level of *ompC* under kanamycin exposure [57]. Moreover, a low level of *ompC* correlates with high kanamycin resistance of the cell.

Our qRT-PCR results based on average Ct values did not reveal any significant transcriptional response of porins to long-term antibiotic exposure for most samples. They said nothing about the variability of *ompF* and *ompC* expression in the *Y. pseudotuberculosis* population or the amount of their gene products in the cell. To answer these questions, we first developed fluorescent reporter strains of *Y. pseudotuberculosis* 488. The GFP reporter systems were constructed by fusion of the *ompF*/*ompC* regulatory regions, including binding sites for transcription and translation factors (5- UTRs and signal sequences for sRNAs, C-terminal YQF amino acids for proteases), to the GFP reporter gene in a low copy number plasmid pACYC184 (Figure S1). Thus, the fluorescence of *Y. pseudotuberculosis* transformants potentially not only provided a high-resolution method for assaying *ompF*/*ompC* promoter activity but also indirectly indicated OmpF and OmpC amounts. We excluded chloramphenicol from the experiments due to the plasmid containing the gene for its resistance. However, we additionally verified that the presence of this antibiotic in the medium did not affect the level of sample fluorescence (Figure S2).

There was little evidence of any change in the mean fluorescence intensity, and hence OmpF and OmpC expression, after prolonged exposure to tetracycline and kanamycin (Figure 5), while qRT-PCR revealed the *ompF* upregulation in the presence of tetracycline at 27 ◦C (Figure 2). This difference indicates an additional post-transcriptional regulation of OmpF.

**Figure 5.** GFP fluorescence intensity in *Y. pseudotuberculosis* 488 transformed with *ompF*/*ompC* promoter-fused GFP reporter constructs following 16 h antibiotic treatment monitored by flow cytometry. (**A**) 27 ◦C incubation temperature; (**B**) 37 ◦C incubation temperature. All results are expressed as mean ± standard deviation between three experimental trials. An asterisk indicates *p*-value < 0.05 vs. respective control. Significance was calculated using one-way ANOVA. The histograms from a single representative experiment are shown. Control—without antibiotic treatment, Km—kanamycin, Tet—tetracycline, and Cb—carbenicillin.

After treatment with sub-MIC of carbenicillin at 27 ◦C, there was a 2.3-fold decrease in the *ompF* level and, on the contrary, a 2.5–3-fold increase in the *ompC* level. This is consistent with the previous observation that adaptation to β-lactams can involve increased expression of OmpC relative to OmpF since an increase in OmpC proportion hinders their penetration [37]. It is worth noting that carbenicillin inhibited the expression of *ompF* only at 27 ◦C. We observed that this transcript was prevailing compared to *ompC* in the untreated samples at both incubation temperatures; however, at 37 ◦C, its level decreased by more than 18 times (Figure S3). This initially low amount of *ompF* transcripts and additional nutrient deficiencies due to the stationary growth phase could prevent a further decrease in its expression at 37 ◦C under carbenicillin stress.

The distribution of GFP fluorescence within the sample cultures under kanamycin and tetracycline exposure was comparable to the untreated controls. Only carbenicillin presence led to the appearance of two *Y. pseudotuberculosis* subpopulations with different fluorescence intensities, most notably at 27 ◦C (Figure 6). F1 subpopulation responded to carbenicillin stress by 1.9-fold inhibition of *ompF* expression, while the cells from F2 increased it by 1.8 times. *ompC* was upregulated under carbenicillin exposure in both subpopulations; however, in C2, its level was significantly higher. Moreover, most of the cells formed subpopulations with a lower level of both *ompF* and *ompC* expression (71% and 78% of the total cell amount, respectively; Figure 6C).

We thus propose that modification of the bacterial envelope by reduced *ompF*/increased *ompC* expression confers adaptive fitness advantages to carbenicillin in a fraction of the bacterial population.

Previously it was shown that heterogenous responses improve the survival of bacterial cultures under the effect of stress factors [58]. Moreover, gene expression variability, including porin genes, is considered to underly adaptive resistance in phenotypically heterogeneous microbial populations [57,59].

**Figure 6.** Fluorescence heterogeneity of *Y. pseudotuberculosis* 488 transformed with *ompF/ompC* promoter-fused GFP reporter constructs following treatment with a sublethal dose of carbenicillin monitored by flow cytometry. (**A**) Fluorescence histograms; (**B**) fluorescent intensity of subpopulations; (**C**) the percentage of cells in each subpopulation. All results are expressed as mean ± standard deviation between three experimental trials. An asterisk indicates *p*-value < 0.05 between groups. Significance was calculated using one-way ANOVA. The histograms from a single representative experiment are shown. F1, F2, and C1, C2—subpopulations of *Y. pseudotuberculosis* 488 with low/high GFP fluorescence.

## **3. Materials and Methods**

#### *3.1. Bacterial Strain and Growth Conditions*

*Y. pseudotuberculosis* 488 strain (= 117 strain, O:1b serotype) isolated from a patient with FESLF in the Russian Far East region was used in the present study.

To study short-term porin response to antibiotic stress, bacterial cells were grown in LB broth overnight at 27 ◦C with shaking. The next day, overnight cultures were diluted at a ratio of 1:50 into fresh LB medium and incubated with agitation at 27 ◦C (optimal growth temperature) and 37 ◦C (host temperature) to mid-exponential growth phase (OD600 = 0.4–0.6). Then, the antibiotics kanamycin (Sigma-Aldrich, St. Louis, MA, USA), tetracycline (Sigma-Aldrich, St. Louis, MA, USA), chloramphenicol (Sigma-Aldrich, St. Louis, MA, USA), and carbenicillin (Sigma-Aldrich, St. Louis, MA, USA) were added to the media in sublethal concentrations (0.5 MIC). After a 1 h antibiotic treatment (180 rpm, 27 ◦C and 37 ◦C), the cells were mixed with RNAprotect Bacteria Reagent (Qiagen, Hilden, Germany) and centrifuged.

To study porin response to long-term antibiotic exposure, overnight cultures were diluted and inoculated into LB broth media to obtain a final concentration of 5 × 105 CFU/mL. Bacteria were incubated with sublethal concentrations of antibiotics (kanamycin, tetracycline, chloramphenicol, and carbenicillin) for 16 h at 27 ◦C and 37 ◦C. RNAprotect Bacteria Reagent (Qiagen, Hilden, Germany) was added to the culture samples.

Three independent bacterial cultures for each test or control condition were prepared as biological replicates for RNA isolation.


*2.4. Classification*

The classification stage employs a multiclass SVM. The basic idea of SVM is to find a hyperplane that separates D-dimensional data into its two classes [38]. SVM is a discriminative model for classification that principally depends on two basic assumptions. First, complex classification problems can be classified through simple linear discriminative functions by transforming data into a high-dimension space. Second, the training samples for SVMs consist only of those data points that lie close to the decision surface, with the supposition that they provide the most relevant information for classification [39]. SVMs were originally proposed as binary classifiers. However, in real scenarios, data is to be classified into multiple classes. This is done by using multiclass SVM. Either a one-against-one (OAO) or one-against-all (OAA) approach can be used [40]. For acoustic scene classification setup proposed in this work, the joint feature vector extracted from previous stage is used to train the multiclass SVM OAO classifier.


**1. Introduction**

354 Remote Sensing – Applications

Mumby, P.J, E.P. Green, A.J. Edwards & C.D. Clark (1997). Coral reef habitat mapping: how

Mumby, P.J. & Harborne, A.R. (1999). Development of a systematic classification scheme of

Mumby, P.J., A.R. Harborne, J.D.Hedley, K. Zychualuk & P. Blackwell (2006). The spatial

Newman, C., E. LeDrew & A. Lim. (2006). Mapping of coral reef for management of marine

Pahlevan, N.; Valadanzouj, M. J. & Alimohamadi, A. (2006). A quantitative comparison to

Philpot, W. D. (1989). Bathymetric mapping with passive multispectral imagery. *Applied*

Primack, R., R. Rozzi, P. Feinsiger, R. Dirzo y F. Massardo. (1998). *Fundamentos de*

Richter, R. (1996). Atmospheric correction of satellite data with haze removal including a

Richter, R. (1998). Correction of satellite imagery over mountainous terrain. *Applied Optics*,

Rivera, J. A. ; Prada, M. C. ; Arsenault, J.L. ; Moody, G. & Benoit, N. (2006). Detecting fish

Secretaría de Marina, Dirección General adjunta de Oceanografía, Hidrografía y

Xin, Y., J. Li y Q. Cheng (2007). Automatic Generation of Remote Sensing Image Mosaics for

Diadema antillarum. *Ecological Modelling*, Vol. 196, pp. 131-148

*Pixels to processes",* Enschede, Netherlands, 8-11 May 2006

*conservación biológica*. Fondo de cultura económica. México.

pp. 88-104. NOAA Professional Paper NMFS, 5 Seattle, WA.

http://www.serc.si.edu/labs/phytoplankton/primer/hydrops.aspx

SERC (September 2011). Phytoplankton lab, 10.09.2011 Available from

72106-2, Springer Berlin Heidelberg New York.

Springer, ISBN 978-1-4020-3968-3, Netherlands

*Optics*, Vol. 28, No. 8 (April 1989), pp. 1569-1578

Vol. 37, No. 16, (June 1998), pp. 4004-4015

7504

326

ISSN 00983004

Scale 1: 100000

1997), pp. 193-202. ISSN: 0025-3162

underwater light fields. *Applied Optics*, Vol. 32, No. 36, (December 1993) pp. 7484-

much detail can remote sensing provide? *Marine Biology*, Vol. 130, No.4 (December

marine habitats to facilitate regional management of Caribbean coral reefs. *Biological Conservation*, Vol. 88, No. 2, (May 1999), pp155-163, ISSN 0006-3207 Mumby, P.J., W. Skirving, A.E. Strong, J.T. Hardy, E.F. LeDrew, E.J. Hochberg, R.P. Stumpf

& L.T. David. (2004). Remote sensing of coral reefs and their physical environment. *Marine Pollution Bulletin*, Vol.48, No. 3-4, (February 2004), pp. 219-228, ISSN 0025-

ecology of Caribbean coral reefs: revisiting the catastrophic die-off of the urchin

protected areas in developing nations using remote sensing. In: *Remote sensing of aquatic coastal ecosystem processes*, Richardson, L. y LeDrew, E., pp. 251-278,

water column correction techniques for benthic mapping using high spatial resolution data. *ISPRS Commission VII Mid-term Symposium "Remote Sensing: From*

haze/clear transition region. *Computers & Geosciences*, Vol. 22, No. 6, pp. 675-681,

aggregations from reef habitats mapped with high resolution side scan sonar imager, In: *Emerging technologies for reef fisheries research and management,* NOAA,

Meteorología (SEMAR) (2008). *Carta náutica SM-932 Majahual to Banco Chinchorro,*

Mapping Large Natural Hazards Areas, In: *Geomatics solutions for disaster management*, Li, J., Zlatanova, S. y Fabbri, A. 451 pp. Springer, ISBN 13 978-3-540-

#### **1.1 Snow reserves and remote sensing**

Snow reserves in mountainous basins are important and reliable water resources in Iran. Identification of their quality is necessary because of an increasing value of freshwater and utilization of water recourses. About 60 percent of surface water and 57 percent of ground water sources in Iran flows in snowy regions (Rayegani, 2005). The water produced from snowmelt process provides soil water, ground water reserves and water in lakes and rivers. Since snow cover is one of the most important sources of provided water, an accurate prediction and timing of snow runoff is necessary for the efficient management and decision- making in water supply.

The science of snow hydrology, compared to other branches of hydrology science, has a relatively shorter history due to difficulties accompanied with snow measurement. The correct analysis of snow issues needs a set of observations and statistics in snow-gauging. Currently, however, there are no regular and comprehensive snow measurement procedures in most parts of Iran. Measurements are only limited to those snowy basins recharging important dams; even these measurements are carried out in scattered points rather than an entire dam catchment area.

The measurement range of these stations is limited to 2000-3000 m asl heights. Thus, in mountainous Iran, current distribution of stations would not seem to be adequate. In such conditions, study of snow reserves and identification of snow melting trend in most basins would be accompanied with limitations. Consequently, measuring snow cover using ground methods will be difficult and costly. Remote sensing technology has many applications in various environmental and earth resources studies including ice and snow research. These applications have been increased recently as a result of unique technical advantages such as multi-temporal imagery acquired in various wavelengths, extent of spatial coverage, and improvement of computer hardwares for interpretation and extraction of information. Regarding snow research, remote sensing technology can provide

<sup>\*</sup> Corresponding Author

Predictability of Water Sources Using Snow

the above-mentioned maps.

or SNOTEL had the highest accuracy.

algorithm will be removed (Taghvakish, 2005; Adhami, 2005).

well, they will reduce the accuracy of snow map.

2. This method is a completely automatic algorithm. 3. This algorithm is applicable for all regions in the globe. 4. This method is simple, accurate and easy to understand.

Maps Extracted from the Modis Imagery in Central Alborz, Iran 357

Snow Map algorithm uses normalized subtractive index (NDSI). Lee et al. (2001) compared MODIS snow maps created with NDSI index with maps prepared by National Operational Hydrologic Remote Sensing Center (NOHRSC, prepared automatically by GOES and NOAA images) in upper region of Rio Grande reservoir. In NOHRSC, the teta algorithm is used. In this algorithm, two classified images are subtracted to identify snow surface. In teta algorithm, two separate threshold limit is introduced for each image. Lee et al. (2001) concluded that both images are affected by cloudy condition and the main error is cloud coverage. They also mentioned that maps produced from MODIS were more accurate than

Ault et al. (2006) concluded that MOD10-L2 snow surface product, MODIS sensor, in clear sky condition had the highest accuracy. They showed that the highest error was associated with those conditions that snow depth was lower than 1 cm; thus the higher was the snow depth, the higher was the accuracy. Hall et al. (2000) also showed a low accuracy of lowmass and patchy snows in New England. Klein and Barnett (2003) carried out a snow cover study using MODIS in Rio Grande reservoir during the 2000-2001 period and compared their results to the ground-measuring methods such as snowpack telemetry (SNOTEL) and NOHRSC models. They ultimately concluded that the highest error associated with maps prepared by MODIS was related to the beginning and end of snowfall period. They showed that when the surface was completely covered by snow with no mixed cloud, ground survey

It can be mentioned that MODIS sensor and NDSI index are appropriate in snow map preparation, although cloud coverage and classification are regarded as constraints (Klein and Barnett, 2003; Zhou et al., 2005). In fact, in spite of various advantages, Snow Map algorithm has some limitations due to inseparability of snow cover from cloud and similarity of cloud behaviours to snow cover. This algorithm cannot completely distinguish clouds from snow (of course, this problem is relatively removable by using Cloud Mask algorithm). Also, this algorithm cannot detect coastal terrains which are similar to snow from viewpoint of whiteness and brightness. However, temperature can act as factor to discriminate snow from these terrains using MODIS bands 31 and 32. Since Cloud Mask and thermal Mask are used before applying algorithm, some error sources in snow map

In snow-gauging using satellite imagery, the existence of cloud is problematic due to the following reasons (Riggs and Hall, 2002): first, clouds conceal Earth information; second, clouds create shades on area and change reflectance. Indeed, if clouds cannot be detected

Clouds and snow have generally similar spectral reflectional properties in range of visible and infra-red spectra, so thermal properties is not enough for distinguishing them from each other as clouds may be cooler or hotter than snow (Singh and Singh, 2001). In order to detect clouds, a procedure called Cloud Mask algorithm is being used. Akerman et al. (1998) introduced MOD35 Cloud Mask algorithm. MOD35 algorithm is based on obstruction of Earth surface affected by cloud or dust particles that identifies water body, land and atmosphere (Strabala,

continuous information layers with higher accuracy and lower cost compared to the ground survey, so it can fill the information gaps in snow hydrological statistics. However, using ground data can increase the efficiency of remotely-sensed measurement of snow-gauging. Satellites are appropriate tools for gauging snow coverage, because of high reflection of snow that creates proper contrast to most of natural surfaces (with the exception of clouds). Therefore, using satellite imagery and GIS modeling one can produce snow-cover maps, assess the changes in snow cover area with various time series, discriminate snow from other features, and model it in a catchment area. These simplify decision-making process for engineers and hydrology managers.

One of the important issues in remotely-sensed snow-gauging is the selection of sensor. Some of optical sensors that have ever used in snow-measuring include sensors mounted to satellites namely TIROS-1 (1960), ESSA\_3, NOAA (1996), LANDSAT (MSS and ETM), and MODIS (2000). Since each sensor has unique properties, a sensor with appropriate spectral, temporal and spatial resolution for snow-gauging must be selected. Since snow is a phenomenon with noticeable surface changes over time, it is necessary to select a sensor that produces proper multi-temporal series. Snow-gauging is done in vast areas, and snow surface is generally even; therefore, MODIS is an appropriate imagery for this purpose. From the view point of spectral resolution, MODIS is one of the best optical sensors for studying snow and discrimination of snow from phenomena such as cloud which has similar spectral reflectance.

One of the purposes of designing of MODIS is a global identification of various types of clouds; hence, several bands have been considered for it to identify various types of cloud cover, optical thickness, effective radius and thermal phase (King et al., 2004).

NASA (National Aeronautics and space administration) launched TERRA satellite to space on December 18th 1999, and MODIS as one of the five sensors mounted on TERRA transferred the first information to Earth on February 24th 2000. MODIS has 36 various bands in visible, infra-red and thermal parts of electromagnetic spectrum including 2 visible bands with 250 m resolution, 5 infra-red bands with 500 m resolution, and 29 thermal bands with 1000 m resolution (Hall et al., 2000).

#### **1.2 A review on remotely-sensed snow measurement**

Various methods have been used to estimate snow surface such as classification methods, threshold limit, decision-based methods, etc. One of the most applicable algorithms used to estimate snow surface is MODIS snow map algorithm. It was introduced in 1998 as a decision-based algorithm which uses group tests of threshold limit for detection of snow. This algorithm has very small volume from the calculation viewpoint and simple from the conceptual viewpoint, thus user can track how product has been created. In addition, this algorithm has an appropriate efficiency with global application (Hall et al., 1998).

Totally the properties of this algorithm include:

1. The precision of this method for various types of snow-covered surfaces for identification of snow surface is higher than other methods such as supervised classification, unsupervised classification and sub-pixel methods provided that atmospheric correction is considered (Dadashi Khaneghah, 2008).

2. This method is a completely automatic algorithm.

356 Remote Sensing – Applications

continuous information layers with higher accuracy and lower cost compared to the ground survey, so it can fill the information gaps in snow hydrological statistics. However, using ground data can increase the efficiency of remotely-sensed measurement of snow-gauging. Satellites are appropriate tools for gauging snow coverage, because of high reflection of snow that creates proper contrast to most of natural surfaces (with the exception of clouds). Therefore, using satellite imagery and GIS modeling one can produce snow-cover maps, assess the changes in snow cover area with various time series, discriminate snow from other features, and model it in a catchment area. These simplify decision-making process for

One of the important issues in remotely-sensed snow-gauging is the selection of sensor. Some of optical sensors that have ever used in snow-measuring include sensors mounted to satellites namely TIROS-1 (1960), ESSA\_3, NOAA (1996), LANDSAT (MSS and ETM), and MODIS (2000). Since each sensor has unique properties, a sensor with appropriate spectral, temporal and spatial resolution for snow-gauging must be selected. Since snow is a phenomenon with noticeable surface changes over time, it is necessary to select a sensor that produces proper multi-temporal series. Snow-gauging is done in vast areas, and snow surface is generally even; therefore, MODIS is an appropriate imagery for this purpose. From the view point of spectral resolution, MODIS is one of the best optical sensors for studying snow and discrimination of snow from phenomena such as cloud which has

One of the purposes of designing of MODIS is a global identification of various types of clouds; hence, several bands have been considered for it to identify various types of cloud

NASA (National Aeronautics and space administration) launched TERRA satellite to space on December 18th 1999, and MODIS as one of the five sensors mounted on TERRA transferred the first information to Earth on February 24th 2000. MODIS has 36 various bands in visible, infra-red and thermal parts of electromagnetic spectrum including 2 visible bands with 250 m resolution, 5 infra-red bands with 500 m resolution, and 29 thermal bands

Various methods have been used to estimate snow surface such as classification methods, threshold limit, decision-based methods, etc. One of the most applicable algorithms used to estimate snow surface is MODIS snow map algorithm. It was introduced in 1998 as a decision-based algorithm which uses group tests of threshold limit for detection of snow. This algorithm has very small volume from the calculation viewpoint and simple from the conceptual viewpoint, thus user can track how product has been created. In addition, this

1. The precision of this method for various types of snow-covered surfaces for identification of snow surface is higher than other methods such as supervised classification, unsupervised classification and sub-pixel methods provided that

algorithm has an appropriate efficiency with global application (Hall et al., 1998).

atmospheric correction is considered (Dadashi Khaneghah, 2008).

cover, optical thickness, effective radius and thermal phase (King et al., 2004).

engineers and hydrology managers.

similar spectral reflectance.

with 1000 m resolution (Hall et al., 2000).

Totally the properties of this algorithm include:

**1.2 A review on remotely-sensed snow measurement**

Snow Map algorithm uses normalized subtractive index (NDSI). Lee et al. (2001) compared MODIS snow maps created with NDSI index with maps prepared by National Operational Hydrologic Remote Sensing Center (NOHRSC, prepared automatically by GOES and NOAA images) in upper region of Rio Grande reservoir. In NOHRSC, the teta algorithm is used. In this algorithm, two classified images are subtracted to identify snow surface. In teta algorithm, two separate threshold limit is introduced for each image. Lee et al. (2001) concluded that both images are affected by cloudy condition and the main error is cloud coverage. They also mentioned that maps produced from MODIS were more accurate than the above-mentioned maps.

Ault et al. (2006) concluded that MOD10-L2 snow surface product, MODIS sensor, in clear sky condition had the highest accuracy. They showed that the highest error was associated with those conditions that snow depth was lower than 1 cm; thus the higher was the snow depth, the higher was the accuracy. Hall et al. (2000) also showed a low accuracy of lowmass and patchy snows in New England. Klein and Barnett (2003) carried out a snow cover study using MODIS in Rio Grande reservoir during the 2000-2001 period and compared their results to the ground-measuring methods such as snowpack telemetry (SNOTEL) and NOHRSC models. They ultimately concluded that the highest error associated with maps prepared by MODIS was related to the beginning and end of snowfall period. They showed that when the surface was completely covered by snow with no mixed cloud, ground survey or SNOTEL had the highest accuracy.

It can be mentioned that MODIS sensor and NDSI index are appropriate in snow map preparation, although cloud coverage and classification are regarded as constraints (Klein and Barnett, 2003; Zhou et al., 2005). In fact, in spite of various advantages, Snow Map algorithm has some limitations due to inseparability of snow cover from cloud and similarity of cloud behaviours to snow cover. This algorithm cannot completely distinguish clouds from snow (of course, this problem is relatively removable by using Cloud Mask algorithm). Also, this algorithm cannot detect coastal terrains which are similar to snow from viewpoint of whiteness and brightness. However, temperature can act as factor to discriminate snow from these terrains using MODIS bands 31 and 32. Since Cloud Mask and thermal Mask are used before applying algorithm, some error sources in snow map algorithm will be removed (Taghvakish, 2005; Adhami, 2005).

In snow-gauging using satellite imagery, the existence of cloud is problematic due to the following reasons (Riggs and Hall, 2002): first, clouds conceal Earth information; second, clouds create shades on area and change reflectance. Indeed, if clouds cannot be detected well, they will reduce the accuracy of snow map.

Clouds and snow have generally similar spectral reflectional properties in range of visible and infra-red spectra, so thermal properties is not enough for distinguishing them from each other as clouds may be cooler or hotter than snow (Singh and Singh, 2001). In order to detect clouds, a procedure called Cloud Mask algorithm is being used. Akerman et al. (1998) introduced MOD35 Cloud Mask algorithm. MOD35 algorithm is based on obstruction of Earth surface affected by cloud or dust particles that identifies water body, land and atmosphere (Strabala,

Predictability of Water Sources Using Snow

coordinates system, PCI Geomatica software was used.

**2.2.2 Digital Elevation Model**

**2.2.3 Ground stations data**

Maps Extracted from the Modis Imagery in Central Alborz, Iran 359

implemented in way that atmospheric diffusion and reflection were minimum. Since data with higher wavelength are being less influenced by aerosols, suspended particles and non-

Digital Elevation Model (DEM) obtained from SRTM Shuttle was used. These data that was in format of GeoTiff had Lambert image system. For transformation of this data to UTM

Snow-gauge station data were obtained from Water Organization- Department of Surface Water. Snow-gauge stations of Central Alborz are located in five basins namely Lattian, Lar, Taleghan, Karaj and Golpayegan basins. In some cases, ground data survey time was not consistent with the time of image acquisition. For solving this problem, in those dates that no ground statistics were available, previous and next day's information were interpolated. Of course interpolation was carried out for those stations where sampling time was close to image acquisition date, and the station had snow cover during the period (February and March).

It is necessary to mention that snow depth data has been used to examine the presence or absence of snow cover. Those gauging stations located far from human interfering features (e.g. buildings) were selected and their snow depth measured. The snow surface was defined an area where surface is regular and even, with the minimum wind blowing effect to increase measurement accuracy (Pfister and Schneebli, 1999). The output spatial resolution of snow map algorithm and Cloud Mask Algorithm is one kilometre. Around each station, up to two kilometres was regarded to include 9 to 13 pixels covering snow. The most repeated pixel shows snowy or not snowy condition. Finally, snow surface obtained from Snow Map Algorithm with and without Liberal Cloud Mask was compared to ground data. Figure 1 illustrates distribution pattern of snow-gauge stations within five basins in Central Alborz.

Among sampled stations, 18 stations with above-mentioned conditions were selected.

Fig. 1. The position Of snow gauge stations in Alborz-e-Markazi

selective diffusion phenomenon, thermal bands of MOD02 imagery were used.

2003). In this process, based on land type, geographical position and available data, Cloud Mask algorithm uses 14 bands amongst 36 bands of MODIS to test 18 spectral and spatial features (Hall and Riggs, 2002). However, this procedure was modified by Hall and Riggs (2002) who presented a new version of Cloud Mask algorithm (Liberal). This algorithm can analyse the pixels located under thin and transparent clouds (Zhou et al., 2005; Ault et. al., 2006). This procedure identifies the darkness and if it faces to such darkness, it means that sun angle is higher than 85º. This algorithm is called Liberal Cloud Mask algorithm. In fact, Liberal Cloud Mask algorithm functions as subset of spectral tests of old Cloud Mask algorithm (MOD35) and uses 7 bands of MODIS and set 4 criteria (Hall and Riggs, 2002, 2004).


**6. Conclusion**

In this study, the accumulated data led to the conclusion that although cognitive decline and brain mitochondrial dysfunction were observed in obese-insulin-resistant rats, AD was not developed during the pre-diabetic state. In addition, markers indicating the presence of AD were observed in T2DM subjects. Treatment with antidiabetic drugs, traditional medicine, FGF21, and mitochondrial fission inhibitors effectively improved brain mitochondrial function and cognitive function in rats with induced MetS.


**3.2. Simulation results**

In this section results of experiments are presented. The purpose of simulation was to determine the optimal number of radio interfaces for different WMN topologies, different number of data flows and different number of nodes to achieve the network capacity increasing expressed in enhancement of QoS parameters.

Channel Assignment Schemes Optimization for Multi-Interface Wireless Mesh Networks Based on Link Load 85

for higher number of radio interfaces. From results it may be concluded that optimal number of radio interfaces which guarantee the maximum allowable average delay 150 ms (ITU-T, 2003) for both network topologies is five, because more than five interfaces improved value of end-to-end delay only slightly, but the complexity of node is increased considera-

Figure 6 shows the average values of network throughput for various numbers of radio interfaces and two different network topologies. The lowest value of average throughput was achieved in grid WMN, where nodes have used for transmission only one radio interface. In this case, the value of average throughput was 504.28 kbps. In the case where WMN with random topology and one radio interface was used, the lowest value of average throughput (739.3 kbps) was achieved. The highest value of throughput (2019.9 kbps) reaches the grid WMN with seven radio interfaces. The best value of average throughput in random WMN topology (1964.2 kbps) was achieved by WMN with seven radio interfaces. Again, the optimal number of interfaces for both network topologies was chosen as five.

**Average Throughput**

Random Grid

**Figure 6.** Average values of throughput for various radio interfaces and different network topologies

where grid topology achieved 9.8% of packet loss and 7.6% for random topology.

As we can see from Fig.7, the highest value of packet loss (75.1%) was reached in grid WMN with one radio interface. The lowest value of packet loss was achieved in WMN with seven radio interfaces. This value was 3.5% for the random topology and 2.5% for grid topology. As in the previous case, we can conclude the optimal number of radio interfaces as five,

12345678 **Number odf Radio Interfaces**

Figure 8 shows the average values of time jitter for different types of topologies and various number of radio interfaces. From results it is obvious that the highest value of average jitter was reached in the network with one radio interface. For the random topology this value was 0.7 sec and for grid topology it was 0.8 sec. On the other hand the lowest values of average jitter were achieved in grid WMN with seven interfaces (0.3 sec) and in random WMN with six interfaces (0.05 sec). As an optimal number of radio interfaces, the number of six was selected with average jitter value 0.11 sec for random topology and 0.14 sec for grid

bly.

topology.

000.00E+0

500.00E+3

1.00E+6

1.50E+6

2.00E+6

2.50E+6

**Throughbut (bit/s)**

We chose four QoS parameters for simulation evaluation:


**5. Workflow**

The following Figure 4 shows the principle work and data flow of the proposed architecture. The ABMS relies on generated data stored in a database provided by synthetic data, Monte Carlo simulation, and data mining (DM) of sensor data from mobile crowd sensing (MCWS). The MCWS is performed by mobile chat bot agents, discussed in Section 6.3, which establishes the connection between virtual and real worlds. But the data flow is bidirectional, and agents can carry data generated by the simulation to mobile devices and users in the real world, e.g., by posting messages in chat blogs or other social media.

**Figure 4.** Principle work and data flow integrating agent-based crowd sensing in agent-based modelling and simulation.

Simulation data is stored in the database, too. The simulation data, e.g., monitoring data or artificial sensor data, can be analysed numerically or statistically by DM. Sensor data can be collected off-line (classical surveys) before the simulation or on-line during the simulation creating incremental simulation runs.


**3. Long-haul tourists**

arrivals in 2015 increased by 4.6% (1180 million worldwide), an increase of 52 million over the previous year. It was the sixth consecutive year of above-average growth in international tourism following the 2009 global crisis [32]. Asia, for example, recorded close to 6% growth in international tourist arrivals. Tourist arrivals have seen growth from 25 million in 1950 to 1.087 billion in 2013 and expected to reach over 1.8 billion in the year of 2030 [32]. According to UNTWO [32], study of motivation of tourists to visit another country, rising economies have become more critical through the fact that their market share increased from 30% in 1980 to 47% in 2013 and is expected to reach 57% by 2013. Nevertheless, based on previous researches on travel motivations [19, 25, 36] in Malaysia, the study elaborated on travel motivations of overall international tourists rather than focused on long-haul

International tourist arrival markets, particularly in the Asian and Southeast Asian regions, are more focused on short-haul tourists than long-haul tourists. The tourism markets in Thailand, for example, received short-haul tourists rather than long-haul tourists. Based on the statistics issued in 2009 [8], Thailand received 21,345,482 short-haul tourists compared to 8,160,880 long-haul tourists (statistics based on the arrival of 25 major markets to Thailand). Likewise, Singapore received short-haul tourists as much as 9,645,346 compared to 3,711,783 long-haul tourists (statistics based on the arrival of 15 major markets to Singapore) [28]. Malaysia also faced the same issue with the arrival of 23,393,445 shorthaul tourists compared to 1,987,545 long-haul tourists (statistics based on the arrival of 15 major markets to Malaysia) [30]. The Penang Institute or SERI [26] said the major tourist markets to Penang in 2014 were from short-haul markets from Indonesia (41.1%) and Singapore (20.4%), compared to long-haul countries, for example, Japan (4.5%), the United States (3.1%), Australia (2.8%) and others. Meanwhile, for tourist entries in 2015 (January– October), only 19.92% (155,263) were long-haul arrivals in Penang compared to 80.08% (779,245) of short-haul tourists [13]. The study was conducted to identify the motivations that affected the long-haul travel by international tourists to Penang and will indirectly fill

Long-haul or long-distance travel can be described as having a long period of time to travel. The concept of long-haul travel based on Dateline project comprises a minimum of 100 km journey [7]. According to Boerjan [3], long-haul travel by airplane involves more than 5 h or travel distance of more than 3000 miles to destination. Letho et al. [16] defined long-haul travel as 'travel more than four nights or more by plane outside of the international area'. In the past, trips involving multiple destinations, long-haul in nature and airplane use were rarely done [20]. Now, travel is already considered habitual, and for some individuals, travel constitutes part of their lifestyle. Introduction to new technologies has enabled people to travel further and comfortably [14–15] (e.g. new airports, direct flights and others) and with further advanced technological developments, travels to a destination are now more conve-

tourists.

the gap of long-haul travel studies for the state.

58 Mobilities, Tourism and Travel Behavior - Contexts and Boundaries

**2. What is long-haul travel?**

nient than ever [16].

Long-haul and short-haul tourists are different based on holidays, trips, budgets, accommodation and others. Most long-haul tourists travel independently [2] and, according to Tidelwall and Faulkner [29] and Hwang and Fasenmaier [12], engage in multi destination trips and sought to satisfy multiple trip purposes. Pre-holiday phase is very important for the long-haul tourists, who usually spend more time gathering information about the trip [3]. According to McKercher and du Cros [18], long-haul tourists were most likely willing to enter a culturally different atmosphere and tend to travel further. Long-haul tourists view their trips as opportunities for learning and personal development as well as seeing it in the context of a rare, often once-in-a-lifetime occurrence [33–34]. Bao and McKercher [2] stated that there were similarities between long-haul tourists and backpackers, such as ability to travel independently, limited budget, simple accommodation and travel pattern. However, based on Crouch [6], long-haul tourists were less sensitive to price and willing to spend more on travel [3], because they consider quality and product features [17].