Iuri Casciuc

Yuliana Zabolotna

Dragos Horvath

Gilles Marcou

Jürgen Bajorath

Alexandre Varnek
email: varnek@unistra.fr

Virtual Screening with Generative Topographic Maps: How many maps are Required?

Keywords: Generative Topographic Mapping, Virtual Screening, Classification Models

Universal Generative Topographic Maps (GTM) provide 2D representations of chemical space selected for their "polypharmacological competence", i.e. the ability to simultaneously represent meaningful activity and property landscapes, associated with many distinct targets and properties. Several such GTMs can be generated -each based on a different initial descriptor vector, encoding distinct structural features. While their average polypharmacological competence may indeed be equivalent, they may nevertheless significantly diverge with respect to the quality of each property-specific landscape. In this work, we show that distinct universal maps represent complementary and strongly synergistic views of biologically relevant chemical space.

Eight universal GTMs were employed as support for predictive classification landscapes, using more than 600 active/inactive ligand series associated with as many targets from the ChEMBL database (v.23). For nine of these targets, it was possible to extract, from the Directory of Useful Decoys (DUD), truly external sets featuring sufficient "actives" and "decoys" not present in the landscape-defining ChEMBL ligand sets. For each such molecule, projected on every class landscape of particular universal map, a probability of activity was estimated, in analogy to a Virtual Screening (VS) experiment.

Calculation of Cross-Validated (CV) Balanced Accuracy (BA) on landscape-defining

ChEMBL data was unable to predict the success of that landscape in VS. Thus, the universal map with best CV results for a given property should not be prioritized as the implicitly best predictor. For a given map, predictions for many DUD compounds are not trustworthy, according to Applicability Domain (AD) considerations. By contrast, simultaneous application of all universal maps, and rating of the likelihood of activity as the mean returned by all applicable Page 1 of 18 ACS Paragon Plus Environment Journal of Chemical Information and Modeling maps, significantly improved prediction results. Performance measures in consensus VS using multiple maps were always superior or similar to those of the best individual map.

Introduction

Nowadays we are facing a growing problem with Big Data in many areas and chemistry is not an exception. Currently, ensemble of academic, commercial and propriety databases records more than 100 million of compounds [START_REF] Polishchuk | Estimation of the Size of Drug-like Chemical Space Based on GDB-17 Data[END_REF] . An estimation of the drug-like chemical space size gives us around 10 33 virtual compounds [START_REF] Polishchuk | Estimation of the Size of Drug-like Chemical Space Based on GDB-17 Data[END_REF] . Hence, selection of potential drug molecules from vast collections of candidate compounds is a real challenge for medicinal chemists.

Chemical information is intrinsically multidimensional, as it may alternatively focus on, for example, connectivity, electronic cloud densities, shape, or pharmacophore patterns, and each aspect may prove to be very important for understanding chemical properties and biological activities. These various properties can be encoded by specific molecular descriptors, i.e. specific vectors of N numbers derived from chemical structure, thus representing a molecule as a point in N-dimensional descriptor space. In principle and at arbitrarily high N, this conceptual space may contain almost all known information about molecules, which, in theory, should allow researchers to predict any desired properties using already obtained experimental values as a training input. However, it is impossible to handle such amount of information without advanced data mining techniques. Even though a variety of methods exist 23 , the main difficulty is striving for a balance between the accuracy of the results and the computational cost of the required calculations.

One of the techniques that is well suited to reach this balance is Generative Topographic Mapping [START_REF] Bishop | The Generative Topographic Mapping[END_REF] (GTM) -a non-linear mapping method that is widely used as visualization tool for analysis of a multidimensional space. GTM landscapes have already been used as QSAR models 567 , and their predictive performance in Virtual Screening (VS) tends to increase with the size and diversity of the dataset used to "color" the landscape. GTM was successfully used for structure-activity analysis of an anti-viral compound set [START_REF] Klimenko | Chemical Space Mapping and Structure--Activity Analysis of the ChEMBL Antiviral Compound Set[END_REF] and also of an anti-malarial mode of action database [START_REF] Sidorov | AntiMalarial Mode of Action (AMMA) Database: Data Selection, Verification and Chemical Space Analysis[END_REF] . Recently, it has also been successfully applied to visualize large public chemical databases such as PubChem, ChEMBL [START_REF] Kayastha | Privileged Structural Motif Detection and Analysis Using Generative Topographic Maps[END_REF] and FDB- [START_REF] Lin | Mapping of the Available Chemical Space versus the Chemical Universe of Lead-Like Compounds[END_REF] . Sidorov et al. [START_REF] Sidorov | Mappability of Drug-like Page 17 of 18 ACS Paragon Plus Environment Journal of Chemical Information and Modeling Space: Towards a Polypharmacologically Competent Map of Drug-Relevant Compounds[END_REF] applied GTM to create "universal" maps of chemical space, that easily distinguished active and inactive compounds for more than 400 ChEMBL targets, yielding an averaged Balanced Accuracy (BA) higher than 0.6 for all targets, indicating high potential of this method for such applications.

The advantage of universal GTM models over classical QSAR approaches is that the most relevant descriptor space that guarantees polypharmacological competence and preferred operational parameter settings defining the manifold are "learned" only once, at the map construction stage. At this stage, large random collections of relevant (drug-like) compounds are used to span biologically relevant chemical space, serving as a "frame set" for unsupervised GTM manifold fitting, while a large and diverse ensemble of structure-activity sets are employed as "selection sets". Their role is to score the quality of the current manifold for its ability to host predictive landscapes corresponding to each selection set activity. Top manifolds scoring well at this stage are selected as the final "universal" maps, with the expectation that they will also be able to support predictive landscapes for other, distinct properties, beyond those present in the selection set. This expectation was well met by more than 400 structure-activity sets consisting of novel compounds associated with completely unrelated targets and properties by Sidorov et al. [START_REF] Sidorov | Mappability of Drug-like Page 17 of 18 ACS Paragon Plus Environment Journal of Chemical Information and Modeling Space: Towards a Polypharmacologically Competent Map of Drug-Relevant Compounds[END_REF] . Certainly, dedicated models that might be built for a given property could exceed the predictive power of universal GTM-based property landscapes -if sufficient training data are available. By contrast, universal GTM manifolds act like "default", zero-parameter models that can even be employed to explore scarcely studied properties with little experimental data. Therefore, they are both the best strategy to use with incipient, small structure-activity series, and an economic, rapid, fitting-free approach to model building for large and diverse series.

In this work, we assess the predictive performance of eight newly constructed universal GTM models in VS of nine target-specific compound sets extracted from Directory of Useful Decoys (DUD) [START_REF] Huang | Benchmarking Sets for Molecular Docking[END_REF] . These GTMs have been constructed on basis of ChEMBL 14 (v.23) structureactivity data for the respective targets -each based on a different initial descriptor vector, encoding distinct structural features. Their average polypharmacological competence is (roughly) equivalent -they are all members of the top ranked population produced by the evolutionary map building process. Nevertheless, they significantly differ in the quality of each property-specific landscape. We show that distinct universal maps represent complementary and strongly synergistic view of chemical space. The predictive power of any classification landscape built for ChEMBL data can be internally assessed by the Cross-validated Balanced Accuracy (BA CV ) criterion in an "aggressive" three-fold cross-validation experiment repeated five times, with data scrambling. However, the BA CV indices were shown unable to predict the success of that landscape in VS. Thus, it would be an error to prefer the universal map with best CV results for a given property as the implicitly best predictor. For a given map, predictions for many DUD compounds are not trustworthy, according to Applicability Domain (AD) considerations. By contrast, simultaneous application of all universal maps, and rating of the likelihood of activity as the mean returned by all applicable maps, significantly improved prediction results. On the basis of different measure, the performance of consensus maps in VS was consistently better than of individual maps.

Methods

Data

The target-specific compound series extraction protocol by Sidorov [START_REF] Sidorov | Mappability of Drug-like Page 17 of 18 ACS Paragon Plus Environment Journal of Chemical Information and Modeling Space: Towards a Polypharmacologically Competent Map of Drug-Relevant Compounds[END_REF] has been applied to release 23 of the ChEMBL database. A total of 618 datasets containing ligands of different ChEMBL human targets have been extracted. The same structure standardization procedure (vide infra) has been applied to DUD database, followed by removal of molecules that were present in ChEMBL to create orthogonal external datasets. For most of the targets shared by ChEMBL and DUD, this required elimination of all the actives from the DUD series. However, in nine cases the DUD target-specific series contained sufficiently numerous original actives and were used for VS. Table 1 summarizes the composition of selected compound datasets.

Workflow

The following workflow was applied:

1) Standardization of ChEMBL and DUD datasets followed by descriptor generation;

2)

Coloring the manifolds of universal maps by each of nine target-specific class landscapes using ChEMBL subsets;

3) 5-fold CV of predictive landscapes using ChEMBL datasets 4) VS applying these landscapes to the DUD subsets;

For some of these steps a dedicated section is presented below. The descriptors used here were ISIDA descriptors computed by ISIDA Fragmentor 161718 .

Data preparation and descriptors generation

More than 100 different types of descriptors sets were generated. They include sequences, atom pairs, circular fragments and triplet counts of different length, colored by formal charges, pharmacophore features or force field types. These fragmentation schemes were selected for relatively low number of fragments they generate.

Generative Topographic Mapping

Generative Topographic Mapping (GTM) is a non-linear mapping method used for data visualization originally described by Bishop. In GTM 2D latent space (called manifold) is embedded into the descriptor space. The points which are close in the latent space remain neighbors in the data space. The manifold represents a grid of k x k nodes; each node is mapped in the initial descriptor space using the mapping function y(x, W). The mapping function is given as a grid of m x m radial basis functions (RBF). In order to build a GTM-based QSAR model, the weighted average of properties of all molecules associated with any particular node is used to "color" the manifold according to that property. Here, the projected property is activity GTM supports several Applicability Domain (AD) [START_REF] Gaspar | GTM-Based QSAR Models and Their Applicability Domains[END_REF] definitions, but only the density-based AD is applied here. Compounds projected onto a "white zone" of the map (accumulating no responsibilities of "training" compounds used to build the landscape) are out of the AD.

Note, however, that the AD considerations in VS may differ from those in predictive QSAR. In the latter case, compounds outside of the AD should be ignored -no prediction of their property should be attempted. In VS, however, the inability to obtain a trustworthy prediction for out-of-AD compounds practically implies that those compounds will be never selected for synthesis and testing because they were predicted to be inactive. Therefore, external compounds falling within the blank spots of the employed class landscapes were assigned zero probability to be active placing them at the bottom of rankings.

Global manifolds (Universal maps) were derived following the procedure in the cited work [START_REF] Sidorov | Mappability of Drug-like Page 17 of 18 ACS Paragon Plus Environment Journal of Chemical Information and Modeling Space: Towards a Polypharmacologically Competent Map of Drug-Relevant Compounds[END_REF] but employing updated compound datasets. They are based on frame sets of maximal diversity (aimed at spanning the entire drug-like chemical space), employed 236 of the abovementioned 618 compound series for map selection. Like in any global mapping approach, they are not meant to capture the detailed SAR of every target-specific set but allow analysis of several activities at the same time. Note that global activity landscapes are relying on a common manifold, itself derived from a selected descriptor space in order to maximize the mean predictive power of all these landscapes. It is obvious that global manifolds represent a best compromise to describe biological activity in general, based on some 'consensus' descriptor space. Interestingly, several such descriptor spaces were identified, each focusing on different aspects of chemical structures. Eight global (universal) maps based on eight distinct ISIDA fragment descriptor spaces were selected (Table 2). On the average, their mean predictive power over all the 618 considered activity sets is similar, while corresponding predictions for each activity series fluctuate. predict the correct activity class of candidates not used for landscape construction, i.e. both in "internal" cross-validation and "external" VS. However, ROC AUC is a more natural VS evaluation criterion than BA, since the latter requires a formal prediction, active vs inactive, for each external compound. In VS, however, the key element is the relative ranking of candidates -a significant prioritization of the actives with respect to the inactives is sufficient to guarantee VS success. Ranking was performed according to the GTM landscape-predicted probability of each compound to be active. The compounds falling outside the applicability domain were assigned zero probability of activity thus they were placed at the bottom of the ranking list. To complement ROC AUC values, the EF of actives ranked within the 100 top compounds was also monitored. EF for the top 100 ranked molecules was calculated according to the equation below.

ܨܧ 100 = ݏ݁ݒ݅ݐܿܣ 100 100 ⁄

ݏ݁ݒ݅ݐܿܣ ݈ܽݐݐ ܰ ݈ܽݐݐ ⁄
where Actives 100 is the number of true positives in the top 100 compounds, Actives total is the total number of active compounds in the dataset, N total is the total number of compounds in dataset.

However, selection of the top 100 compounds may be considered only if there is a significant gap between the probabilities to be active of the 100 th selected compound and the Page 8 of 18

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling one of the 101 st , not selected candidate. In practice, several candidate compounds will have the same predicted probability to be active (reported with a precision of 0.01) and therefore all those that are equiprobable to the 100 th selected compound would be equally deserving to enter the selection. In order to force selection of a top 100 compounds, a random subset of these equiprobable must be picked in completion of the better ranked candidates. In this a posteriori study, three scenarios are considered to compute the EF:

1. Pessimistic: out of candidates that are equiprobable to the 100 th selected compound, inactives are selected first, and then the remaining places in the pessimistic top 100 are completed by actives.

2. Optimistic: the opposite strategy (actives are filled in first, remaining places taken by inactives)

3. Stochastic pick out of candidates that are equiprobable to the 100 th selected compound.

Scenarios 1 and 2 are deterministic. The values obtained are termed Pessimistic Enrichment Factor (PEF) and Optimistic Enrichment Factor (OEF), respectively. Scenario 3 is not deterministic and repeated random drawing/averaging would be required to converge to expectation values. Yet, it is possible to estimate an average value, termed Interpolated Enrichment Factor (IEF) using the following equation:

ܨܧܫ = ߣ × ܨܧܲ + ሺ1 -ߣሻ × ܨܧܱ ߣ = ݊ ܰ
where IEF -interpolated enrichment factor; OEF -optimistic enrichment factor; PEFpessimistic enrichment factor; λ -the ratio n/N, with N being the size of set including all the candidates that are equiprobable to the 100 th selected compound and n the number of these latter candidates. For instance, if the set including all the four candidates that are equiprobable to the 100 th selected compound contains 102 hits, then N=4and n=2 such that λ=0.5.

Results

Cross-Validation of ChEMBL activity class landscapes

Three-fold CV of the BA was repeated five times for each of the ChEMBL series. For the 236 "selection" series, this was part of GTM manifold scoring process, where the fitness score reflects the mean of each BA CV value. For the eight selected manifolds, the same CV procedure was applied to the remaining 618-236 "external" series, thus obtaining the complete matrix of the predictive power of every map for each of the 618. Unsurprisingly, not every property is equally well predicted by each map, albeit the average BA CV value may not differ much from map to map. Each map was examined in order to identify the number of targets for which it is able to solve the active/inactive classification problem at BA CV above a given threshold. Figure 3 shows that for 617 of 618 targets, BA CV scores of 0.6 or better are achieved by at least one of the maps. The exception (ChEMBL5678) represents a set with too few compounds. Note that maps are ranked according to their original fitness score (mean BA CV scores over the 236 selection SAR series) and it can be seen from Figure 3 that the first map is strongly predictive (BA CV >0.75) for 418 distinct series. Note that part of these 418 are selection series but include a significant number of external series nevertheless. It is also noteworthy that every single map is able to provide significantly better-than-random separation of actives and inactives (BA CV >0.6) for virtually all (609/618 -in case of map #1) SAR sets, which fully justifies the label of "universal" maps. However, one single map is not expected to flawlessly model all series -no single descriptor space (fragmentation scheme) on which a map is built could capture all the relevant chemical information that might impact so many different structureactivity relationships. The eight selected maps are highly complementary: series less well explained by one map will work better on another manifold, exploiting specific information from its distinct descriptor space to host a strongly predictive model. Cumulated prediction performance increases with the number of considered maps (Figure 3) which clearly demonstrates the maps complementarity: Seven universal maps based on as many distinct descriptor spaces are sufficient to provide at least one satisfactory result for more than 85% of used targets even at the very stringent BA CV >0.75. Thus, for further analysis, only seven universal maps were used. Next, the question how to identify the best universal map for a particular activity was addressed. It may be expected that the model which shows highest predictive CV performance in target-specific ligand classification would be the best model in VS. In order to test this hypothesis, correlation between landscape performance in CV and VS was evaluated for each of the 63 QSAR models (activity landscapes for nine targets on seven universal maps). Figure 4 compares, for the specific activity landscapes of target CHEMBL260 hosted on each map, the "internal" estimation predictive power (BA CV ) on one hand, and, the observed predictive power in "external" VS of the DUD subset, on the other hand. The Pearson correlation coefficient of BA CV versus BA VS over the seven maps was calculated for all the nine sets; they vary in the range 0.02 -0.63 which means that a map can hardly be chosen on the basis of its CV performance. Unfortunately -but not unexpectedly [START_REF] Golbraikh | Beware of Q2![END_REF] high BA CV is a necessary, but not sufficient guarantee of model success in VS. The success in a predictive challenge depends on the peculiar composition of the test set.

Consensus of Universal Maps.

Given the genuine complementarity of the seven maps, consensus predictions by averaging results these complementary views of chemical space might be a promising strategy.

Here, averaging was applied to the probability of activity from each of the seven landscapes, for each compound from the external test set, excluding, however, landscapes in which the compound was projected into an "empty" zone (Figure 5). In this study, the density-based AD criterion as implemented by default in ISIDA GTM)was applied [START_REF] Gaspar | GTM-Based QSAR Models and Their Applicability Domains[END_REF] . Compounds that fell outside the AD in all existing maps were excluded from the consensus model. Apart from the fact that consensus allows making predictions without choosing a priori one best map, it has another important advantage -data coverage increase (percentage of the compounds that are considered to be in AD). For example, none of the maps of the CHEMBL260 subset provided 100% data coverage achieved by the consensus. Similar observations were made for the remaining eight datasets. Only for two, coverage was less than 100% (ChEMBL4338 -79,8%; ChEMBL4439 -97,5%). Recall that in a VS context, compounds out of AD are not "discarded", but given a probability of zero to be active, which implicitly ranks them at the bottom of the list. Thus, data coverage in this context does not impact on the size of the screened compound set (BA, EF and ROC AUC values are reported with respect to the full In Error! Reference source not found. the results for all nine datasets are shown. The consensus model performed than any individual map on the basis of EF. We note that CHEMBL4338 represents an atypical dataset because about 60% of the compounds fell outside the applicability domain. Analysis of CHEMBL4338 revealed the presence of distinct structural subsets in DUD and ChEMBL, which provided a rationale for low performance on the basis of BA, ROC AUC and EF.. The ChEMBL series used to build the activity landscape mainly contained fused aromatic heterocycles such as hypoxantine, pyrrolopyrimidne, benzimidazole-4,7-quinone (Figure 7). In the DUD series, the majority of compounds that were correctly predicted contained a purine moiety similar to training set molecules. (Figure 7). However, compounds outside of AD were N-phenylsulfonamides (Figure 7) that were not present in the ChEMBL dataset.

Conclusion

A new series of "universal" chemical space maps from datasets in the ChEMBL23 database was built using the GTM dimensionality reduction algorithm and following a previously reported evolutionary procedure to select preferred descriptor spaces and GTM parameter stings. These maps were able to provide better than random separation (BA CV >0.6) of actives and inactives in 609 of 618 ChEMBL sets, irrespective of whether series were used for map selection or not. However, consistently accurate predictions for each activity class could not be achieved be achieved by any individual map. However, these maps, which were each based on a different descriptor space, were highly complementary. For 617 of 618 activity classes, at least one out of the seven top universal maps represented a highly discriminatory activity landscape.

Since there is no correlation between performance in CV and external predictive power of individual activity landscapes, the one possible solution is to use a consensus approach. The, all landscapes with favorable density distributions of VS candidates make positive contributions to the consensus model. The most important advantages of a consensus map are: 1) 100% data coverage in most of the cases; 2) significant increase in EF for the 100 top ranked compounds; 3) high performance of the consensus model compared to individual models on the basis of ROC AUC. Thus, while any single universal map displays moderate predictive power f, the combination of complementary maps results in a strong consensus effect in VS. Seven universal maps were sufficient to generate complementary views of biologically relevant chemical space that resulted in further increased VS performance.

Supporting Information

Activity landscapes for all nine DUD subsets used in VS are provided.

••

  Structures from both databases ChEMBL (version 23) and DUD were standardized accordingly to the procedure implemented on virtual screening server of the Laboratory of /webserv/VSEngine.html) using the ChemAxon Standardizer 15 : Dearomatization and final aromatization according to the "basic" setup of the ChemAxon procedure (heterocycles like pyridone are not aromatized) • Dealkalization • Conversion to canonical SMILES • Salts and mixtures removal Page 5 of 18 ACS Paragon Plus Environment Journal of Chemical Information and Modeling • Neutralization of all species, except nitrogen (IV) Generation of the major tautomer according to ChemAxon After the standardization, 1 540 615 compounds from ChEMBL and 914,379 compounds from DUD remained.

  class membership, resulting into a fuzzy activity landscape. Molecule "responsibilities" are used as weights. Red and blue zones are only populated by active and inactive compounds, respectively; all colors in between correspond to the regions occupied by compounds of both classes in different proportions. White zones represent unpopulated areas.

Figure 1 :

 1 Figure 1: A frame set of compounds is represented in the N-dimensional descriptor space. A flexible 2D manifold,which is a square grid of nodes, is injected into that space and is fitted to the data. The molecules are non-linearly projected onto it, and when the manifold is unbent, a 2D map is obtained. Each node can be colored according to the activities of molecules residing there, producing "activity landscapes", where red zones are populated only by active molecules, blue -by inactive, all colors in between correspond to the regions occupied by compounds of both classes in different proportions. White zones are empty.

Page

  fragments of restricted atom and bonds of a length 1 to 2 atoms labeled by a Force Field.3172 3 IAB-PH-FC-AP-2-4 Sequences of atoms and bonds of a length 2 to 4 atoms labeled by 4245 Model performance was evaluated using BA in 3-fold CV and VS, Receiver Operating Characteristic Area Under Curve (ROC AUC) in VS and Enrichment Factor (EF) in VS. BA has been mainly used during cross-validation. BA serves to assess the ability of landscapes to

Figure 2 :

 2 Figure 2: Heatmap showing the performance of universal maps on 618 selected series. Color-codes: dark blue -BA>0.85, light blue -0.65 < BA ≤ 0.85, orange -0.5 < BA ≤ 0.65 and red BA ≤ 0.5. Between parenthesis is shown the number of target-specific classification problems for which a map scores BA > 0.75.

Page 10 of 18 ACS

 18 Paragon Plus Environment Journal of Chemical Information and Modeling

Figure 3 :

 3 Figure 3: Cumulated performance of universal maps: number of predicted target-specific series vs number of used maps

Figure 4 :

 4 Figure 4: BA values obtained in CV and VS of the CHEMBL260 dataset.

Figure 5 :

 5 Figure 5: Activity landscapes of the CHEMBL260 dataset in seven universal maps.

Page 13 of 18 ACS

 18 Paragon Plus Environment Journal of Chemical Information and Modeling DUD sets, respectively). Data coverage however impacts the reliability of results since increasing data coverage reduces compounds with zero probability of activity.

Figure 6 Figure 6 :

 66 Figure 6 shows that consensus BA values generally exceed the majority of BA scores achieved by individual universal maps. Only universal the map of CHEMBL260 outperformed the consensus model

Figure 7 :

 7 Figure 7: Maximum common substructures of compound subsets active against the purine nucleoside phosphorylase receptor in the CHEMBL4338 dataset and DUD.

Table 1 : Description of target-specific subsets used for model training (ChEMBL) and VS (DUD).

 1

	ChEMBL ID	Target Name	DUD dataset Active Inactive	ChEMBL dataset Active Inactive
	1827	Phosphodiesterase 5A	170	25334	691	1515
	1952	Thymidylate synthase	63	6113	124	455
	251	Adenosine A2a receptor	79	28001	1303	3618
	260	MAP kinase p38 alpha	100	32925	1453	2567
	279	Vascular endothelial growth factor receptor 2	94	22595	2047	4663
	301	Cyclin-dependent kinase 2	189	25675	638	2305
	4282	Serine/threonine-protein kinase AKT	52	14228	725	2619q
	4338	Purine nucleoside phosphorylase	102	6334	100	111
	4439	TGF-beta receptor type I	82	8013	282	385

Table 2 :

 2 Description of eight universal maps, their descriptor types and the descriptor space dimensionality.

Table 3 :

 3 Performance in CV and VS for individual universal maps compared to consensus models.

		Cross-Validation		Virtual Screening			Consensus model
	Target	Best Map number	BA	Best Map number	BA ROC AUC	EF	BA	ROC AUC	EF
	ChEMBL1827	4	0,82	7	0,70	0,73	0,00 0,67	0,74	1.5
	CHEMBL1952	4	0,83	5	0,82	0,85	0,13 0,82	0,86	14.7
	CHEMBL251	2	0,77	3	0,77	0,84	1,56 0,80	0,88	17.8
	CHEMBL260	2	0,75	5	0,71	0,73	0,00 0,64	0,77	11,00
	CHEMBL279	2	0,73	4	0,71	0,78	0,00 0,66	0,82	4.83
	CHEMBL301	3	0,80	5	0,74	0,80	0,60 0,81	0,87	5.47
	CHEMBL4282	5	0,81	3	0,81	0,87	17,39 0,83	0,92	52.18
	CHEMBL4338	5	0,84	3	0,71	0,73	0,00 0,54	0,66	0,00
	CHEMBL4439	5	0,81	5	0,75	0,88	1,97 0,67	0,88	4.94

Acknowledgment

IC thanks the Région Grand Est for a PhD fellowship. Page 16 of 18 ACS Paragon Plus Environment Journal of Chemical Information and Modeling


Marta Glavatskikh

Timur Madzhidov

Dragos Horvath

Ramil Nugmanov

Timur Gimadiev

Daria Malakhova

Gilles Marcou

Alexandre Varnek

Predictive Models for Kinetic Parameters of Cycloaddition Reactions

Keywords: cycloaddition reactions, QSPR, Condensed Graph of Reaction, Generative Topographic Mapping

come   L'archive ouverte pluridisciplinaire

Introduction

Chemical reactions represent rather complex objects for structure-property modeling because their yield, thermodynamic and kinetic parameters depend not only on structure of reactants, but in great extent on experimental conditions. Therefore, an effort should be made to encode the above information in molecular descriptors vector which serves as input in machine-learning methods.

This paper is devoted to development of predictive models for three kinetic properties of cycloaddition reactions: rate constant logarithm (logk), activation energy (Ea) and preexponential factor logarithm (logA). Cycloaddition (CA), a classic reaction in organic chemistry, is an important tool to produce compounds of cyclic architecture, with high regio-and stereoselectivity. In early works, frontier molecular orbitals (FMO) [1][2][3][4][5] calculated by quantum mechanics methods for small congeneric series were used for interpretation of thermodynamic and kinetic properties of cycloaddition. An effort has been made to build a linear correlation with experimentally measured physicochemical parameters of the reagents [6][7] . This, however, significantly limits application of the latter to already studied reactions.

Classical QSPR modeling of chemical reactions was performed on homogeneous series, either varying the reagents under fixed reaction conditions (solvent, temperature), or inversely varying conditions to study a reaction between given reactants [8][9][10][11] . In these studies topological indices [12] , quantum-chemical [8][9][13][14] or mixed [15][16][17] descriptors for reagents were used. For an non-exhaustive overview of the studies carried out in the field of chemical reactivity one can refer on papers by Warr [18] , Baskin et al., [19] or on earlier works [10][11]20] .

In spite of significant progress in the field, few attempts have been made towards a more general QSPR treatment of chemical reactions. Thus, Marcou et al [21] reported classification models predicting experimental conditions (solvent type, catalyst) of the Michael addition.

Predictive models for the rate constant of SN2 [22][23] and E2 [24][25] reactions explicitly accounting both structure of reactants / products and experimental conditions have also been reported in our early publications.

In this work, we employ the approach suggested in [22][23][24] to building predictive models for three kinetic properties of cycloaddition reactions: logk, Ea and logA. The transformation of reactants into products was encoded by a Condensed Graph of Reaction (CGR) [26][27] representing a given reaction as a single pseudomolecule (Fig. 1). CGR-based descriptors thus explicitly characterize the reaction center together with its immediate neighborhood. In our previous works, the CGR-based description has been successfully used for the prediction of rate constant of SN2 and E2 reaction [24][25] , tautomerization equilibrium constants prediction [28] , and the prediction of optimal conditions [27] . Here, the CGR-based descriptors and the reaction conditions descriptors have been concatenated thus forming a descriptors vector for chemical reaction. These vectors were used as an input in model building with Support Vector Regression (SVR) [START_REF] Drucker | Advances in neural information processing systems[END_REF] and Generative Topographic Mapping (GTM) [START_REF] Bishop | [END_REF] machine-learning methods. GTM as a probabilistic extension of self-organizing maps, is an efficient method of visualization, analysis and modeling of chemical properties [31][32] or biological activities [33][34][35] and could be used for both classification and regression tasks. The algorithm converts the N-dimensional data into a 2D map that could be 'colored' by different reference properties. Here, the same 2D map built for the ensemble of cycloaddition reactions was used for the analysis of the logk, Ea and logA parameters.

Computational procedure

The modeling workflow is shown in Figure 2

Data preparation

An initial 2551 cycloaddition reactions with available experimental values of the rate constant (logk), 1356 reactions with activation energies (Ea, in kJ/mol) and 1237 with pre-exponential factor (logA) values were collected manually from the manuscripts of PhD thesis works of Prof. Konovalov's group from Kazan Federal University, published in 1970-1990. Within the 2551 reactions there were 806 unique structural transformations. The data set contained about 85% of Diels-Alder (4+2) cycloaddition, about 8% (3+2) dipolar cyclisation, and 7% (2+2) cycloadditions. The reactions were carried out in C6H5CH3, CH3COOC2H5, C6H5OCH3, C6H5NO2, CH3CN, C2H5OH, C6H5Cl, THF, CH2Cl2, DMSO, CHCl3, 1,4-Dioxane, C5H11OH, C6H5Br, C6H6, CH3COCH3, C6H12 and CH2ClCH2Cl, at temperatures ranging from 273 to 423 K. The most frequently occurred dienes for the Diels-Alder series are: condensed aromatic cycles, cyclopentadienones, cyclopentadienes and benzofurans. Among dienophiles, the prevailing structures are: maleic acid derivatives, cyanoethylenes, Ph-substituted ethylenes and benzoquinones. The data set (reagents and products) underwent an accustomed standardization protocol using ChemAxon's Standardizer Utility: basic aromatization, isotopes removal, NO2-, NO-, N3-, RRSO2-, CN-group transformation [36] .

Cyclization may lead to different stereoisomers, and available data displayed several examples of logk differences associated to the competing processes, all within 0.01 to 1.0 log values. Based on our previous work [24] we considered that such difference is too small to be taken into account, and comparable with the interlaboratory errors affecting logk measurements. Thus, for stereospecific reactions, the logk constants were calculated as a mean of given experimental values. No other duplicates were found in the dataset. The 1849 reactions kept after

Descriptors

The reactions were first rendered as Condensed Graphs of Reaction (CGR), created from the reaction RDF files using the in-house CGR Designer tool and stored in modified SDF format. This was directly processed by the in-house ISIDA Fragmentor software, in order to generate fragment descriptors. The length of monitored fragments varied from 2 to 14 for sequences and from 2 to 6 for atom-centered fragments. The following options were also used at choice: charges on atoms (Formal Charge), accounting for the terminal atoms of a fragment exclusively (Atom Pairs), and exploring all possible paths instead of shortest paths (DoAllWays). An important option regulating the amount of the overall generated CGR fragments is the 'dynamic bond' exclusive inclusion. Toggled on, the option produces the fragments, that contains the bonds forming/breaking while chemical reaction (local fragments) and omits the 'generic' fragments, not assigned to the reaction center. That could be used to generate fragments that describe local environment of the reaction center exclusively. For this project, the CGR fragmentation implied the generation of all possible, local and nonlocal, fragments. Overall, 728 descriptor sets have been generated for the preliminary SVR scanning. Descriptors of the reaction conditions. Descriptors of reaction conditions included solvent descriptors and temperature. The solvent descriptors considered in this study include the values of polarity, polarizability, H-acidity and basicity: Catalan SPP [37] , SA [38] and SB constants [39] , Camlet-Taft α [40] , β [41] , and π* [42] constants, 4 functions of dielectric constant ε (Born

1 B f and Kirkwood 1 2 1 K f functions, 1 1 1 f , 2 1 2 f
), 3 functions of the refractive index nD 20 (denoted as n for the sake of simplicity in the following formulae), 2

1 2 2 1 n n g , 1 2 1 2 2 2 n n g , ) 1 2 )( 1 2 ( ) 1 )( 1 ( 2 2 n n h
). The final element of the condition descriptor features the reciprocal of reaction temperature, given in Kelvin -1 (K).

Building and validation of the models

SVR modeling. SVR models were built and validated using the ε-SVR algorithm implemented in the libSVM package [43] . The modeling was performed using the evolutionary SVR optimizer [44] ,

which can be used to perform both descriptor space selection and optimization of the operational parameters (epsilon, kernel type, cost, gamma) of the SVR method. The procedure, applied to logk as a modeled property generated a total of 6000 models. The descriptor set producing the SVR individual model of maximal robustness (estimated by five-fold cross-validation Q 2 value) was based on atom-centered fragments of the length 2-3. This descriptor space encompassing 688 descriptors (674 of fragments, 14 of the reaction conditions) was used further for SVR and GTM model building and its external validation.

GTM modeling. GTM constructs a 2D-dimensional map reproducing the data distribution of the initial D-dimensional space, defined by the descriptors. The map is framed by K nodes forming a perfect square grid. The nodes correspond to normal probability distributions (NPD) centered on the GTM manifold, a flexible 2D envelope embedded into the D-dimensional data cloud. The assignment (mapping) of the nodes to the manifold points is defined by a mapping function set up with the help of M radial basis functions (RBFs):

where d goes from 1 to D, W is the weight matrix connecting Gaussian RBF and data space points, is the center of the m-th RBF in the 2D space, and is a width of RBF. The nodes framing the GTM map thus correspond the centers of RBFs coordinated on the manifold so that to best approximate the data distribution. The RBFs determines the probabilities of an object to be mapped into a certain node. Moreover, the data object has non-zero probability over all the nodes. Consequently, the object could be characterized by its probabilities, which are called responsibilities , that represent the degree of association between the object n (in the initial D-dimensional data space) and the node k. Details of the GTM algorithm are given in our publications [31,[33][34][35] and references therein. GTM property landscapes. Distribution of the property over the chemical space can be visualized as a GTM property landscape, in which the mean property value of each node of the grid defined in the latent 2D space can be computed as:

where is the property value of the nth molecule, are responsibilities of the molecules to be located in the node , and N is the number of molecules in the data set. If a property landscape in rendered in 2D, the values of of the k-th node determines its color, while the transparency depends on its occupancy (to be more exact, cumulated responsibilities): the more molecules are located near a certain node, the opaquer is the color.

The GTM models were built using the evolutionary optimizer [44] which supports the choice of the operational GTM parameters (the number of RBF kernels, the number of grid points, the width factor of radial basis functions, and the regularization coefficient). The descriptor space, based on atom-centered fragments of the length 2-3, chosen by the preliminary SVR scanning was considered. The genetic optimization procedure included 3000 generations.

The current generated GTM manifold is used as "support" for all the three logk, Ea, logA property landscapes, which serve as predictive models and undergo five-fold cross-validation.

The fitness score of the current GTM manifold (the objective function to be maximized) represents the mean <Q 2 > value of scores achieved by each property landscape. The most robust GTM model, in the above-mentioned sense of maximal <Q 2 > has been further used for the prediction and visualization of logk, Ea and logA.

Validation of the models. The predictive performance of the GTM and SVR models has been estimated by root mean squared error (RMSE) and squared determination coefficient calculated in five-fold cross-validation (Q 2 ) repeated 10 times after the data reshuffling (10x5-CV), or on the external test set (R 2 ):

Here Yexp and Ypred are, respectively, experimental and predicted values of log k, n is the number of data points, while <Y>exp is the mean of experimental values. Applicability Domain of the models. The "Fragment Control" [45] procedure was chosen to delimit the general Applicability Domain [46] (AD) of the approach, which corresponds to the bounding box encompassing the training set in fragment descriptor space. For each of the fragment associated to descriptor vector elements, a candidate molecule must not contain more of these than the richest representative of the training set, or less than the poorest one. Presence of fragments never seen among training compounds therefore counts as AD violation. Since both SVM and GTM models are based on the same fragment descriptor space, their AD are identical.

Different scenarios of logk assessment for the test set reactions

The logk values for the external set of 200 reactions were obtained using three different approaches: direct assessment, Arrhenius-based assessment and temperature-scaled GTM landscapes. The former corresponds to the application SVR or GTM predictive models for logk to the test set reactions.

In the Arrhenius-based assessment the logk values are calculated using the Arrhenius equation (eq. 5). The latter implies the usage of predicted Ea and logA values obtained with the help of related SVR or GTM models. In the Ea and logA modeling, temperature descriptor was excluded.

The rate constant may be very sensitive to temperature. This effect might not be properly accounted for by GTM-based models, where temperature is just one of hundreds of dimensions controlling the manifold geometry. Therefore, instead of using temperature as a descriptor, the temperature-scaled approach has been proposed. It implies construction of a series logk GTM landscapes each corresponding to a specific temperature range.

Scaling of logk measured at temperature T1 to temperature T2 can be performed according to the van't Hoff equation:

. The temperature coefficient was computed as an average of (i = 1, 2… n) values calculated for n series of the reactions, each studied at several temperatures. Overall, reaction rate of 358 reactions in the training set was measured at 2 -6 different temperatures. The smallest temperatures difference in a series was 10K and the largest was 110K. The temperature coefficient varied from 0.08 to 0.59; its average value 0.26 was used in our calculations.

The entire temperature range (273K-423K) has been divided into 15 subranges of 10K each, e.g. from 273K to 283K. For each subrange, the corresponding logk landscape built for all 1649 reactions was constructed, i.e., each logk was recalculated to median temperature of a subrange (e.g., 278K for the range 273K-283K) using van't Hoff's equation. Notice that the logk error related to a temperature deviation of 5K (0.15-0.3 logk units) is far less than related RMSE values (Table 2). Thus, at the validation stage, the program selects a particular logk landscape corresponding to the external reaction temperature.

Thus, fifteen temperature-scaled logk landscapes were calculated from experimental logk values using the Van't Hoff relationship, and related median temperature of a subrange. These landscapes were used for validation of GTM-based models on reactions from the external test set.

Results and discussions

GTM visualization of the training set

The Generative Topographic Map shown on Fig. 4 shows well separated zones populated by the (4+2), (3+2) and (2+2) cycloaddition reactions. The ratio of sizes of these zones correspond to class occurrences in the data set: thus, the major class (4+2) occupies the largest area, flowed by the (3+2) zone, which in turn is slightly bigger than the zone populated by (2+2) reactions. A distinct separation of the cycloaddition classes on the map infers a competency of the chosen descriptors space, able to correctly discriminate the reactions belonging to different types. Another technique of GTM visualization is property landscapes, characterizing the distribution of a certain property over the training set projection. Figure 5 shows GTM property landscapes for logk, Ea and logA, where the gradient from blue to red corresponds to the property variating from low to high value (in the range its own for each property). Simultaneous analysis of all three landscapes helps to understand decomposition of logk on Arrhenius equation parameters logA and Ea (eq. 5). According to collision theory, A is the frequency of collisions in the "correct orientation". Thus, logA values might be related to steric interaction of reactants in a pre-reaction complex. Activation barrier Ea accounts for electronic and steric effects in the transition state. Typically, high reaction rates correspond to large logA and small Ea (zone A on the maps) whereas low logk values result from a combination of small logA and large Ea (zone B). However, in some cases low activation barriers do not compensate the very low logA (zone C), or, on the contrary, large logA is compensated by high activation barriers (zone D). Intuitively, the example C in the Table below is a clear example of sterically crowded transition state (low logA), not compensated by the otherwise favorable electronic interactions between an electrondepleted diene and a neutral dienophile. Example D is, by contrast, featuring a less favorable interaction of electron-depleted diene and electron-depleted dienophile, however not involving any particular steric crowding.

Cross validation of the SVR and GTM models

The cross-validation results (5CV) for the three properties, logk, Ea (kJ/mol) and logA, for SVR and GTM models are shown on Table 2. The SVR modeling involved the construction of three different models for each property, whereas for the GTM, they were assessed by three property landscapes based on one same manifold (see 2.3). One may see that performances of SVR models exceed those of GTM model -unsurprisingly, knowing that the molecular descriptor set best suited for SVM modeling was the one used throughout this work. However, the overall statistics for the GTM is reasonable for all three properties.

External validation of SVR and GTM models

The external test set has been predicted in three different ways (see 2. -(iii) prediction of logk using temperature-scaled GTM landscapes.

The results for both strategies are given in Table 3 and Figure 6. Similarly to cross-validation results, SVM models perform better than GTM ones. The models with direct logk assessment lead to better results than those involving the Arrhenius equation, which might be explained by both smaller training sets for logA and Ea modeling, and the fact that some test set reactions are dissimilar compared to those in the training set. Indeed, the Fragment Control" [45] applicability domain discards 36 out of 200 reactions which leads to significant improvement of models' performances (Table 3). 3.

Reaction rate assessment with temperature-scaled logk landscapes shows a minor improvement over the direct GTM logk prediction (Table 3). This trend was confirmed in logk predictions for the specific subset of 13 reactions measured at two different temperatures (see Fig.

Subset of new transformations.

Since the external test set has been randomly chosen from the initial data set, it includes several reactions present in the training set but associated to different conditions. Therefore, prediction propensity with respect to novel transformations and to different conditions for a known transformation could be assessed separately. Among 200 external test set's reactions, only 57 transformations didn't occur in the training set. Among them, 8 reactions were found out of the AD. Table 11 shows that only direct logk assessment with SVM models leads to reasonable agreement between predicted and experimental values for 49 reactions retained by AD. Indeed, the model performs similarly to cross-validation (R 2 = 0.87, RMSE = 0.59 log units). On the other hand, application of the Arrhenius-based SVR assessment as well as GTM-based models result in poor predictions with R 2 < 0.5. 4.

Conclusion

In this paper we report the first SVR and GTM models of kinetic parameters (reaction rate, preexponential factor in the Arrhenius equation and activation energy) built on large and structurally diverse data set of cycloaddition reactions. A reasonable accuracy of reaction rate assessment was achieved both in cross-validation and on external test set (RMSE = 0.45 and 0.35 log units, respectively). Different reaction classes are well separated on 2-dimensional Generative Topographic Map. The GTM-based property landscapes provide with an interesting insight on reaction chemical space and allow to identify different types of reactions with respect to the interplay between the Arrhenius parameters logA and Ea. The obtained results demonstrate an efficiency of the CGR-based fragment descriptors for chemical structures, and a proper choice of the parameters encoding experimental conditions. The SVR models for the rate constant, activation energy and the pre-exponential coefficient are freely available on our web-server: http://cimm.kpfu.ru/models.

Figure 1 .

 1 Figure 1. Reaction of cyclization and the corresponding CGR of the transformation. Formed bonds are denoted with circle and the broken ones are crossed.

  . The database of 1849 reactions of cycloadditions, for which the logk, Ea and logA values have been measured in different solvents and at different temperatures, has been collected from the literature (see 2.1). The data set has been divided into training and test sets, where the latter has been composed out of 200 reactions picked up randomly from the data set. The remaining 1649 reactions constitute the training set, which has been used to build three individual models, for logk, Ea and logA prediction, respectively. The obtained best SVR and GTM models have been further challenged for the exhaustive external prediction of the logk on the test set of 200 reactions.

Figure 2 .

 2 Figure 2. Workflow for the modeling of the rate constant (logk), activation energy (Ea) and pre-exponential coefficient (log A) of the reactions of cycloaddition.

  cleaning have been characterized by 1849 values of logk, 1356 values of Ea (kJ/mol) and 1237 values of logA, the latter being a strict subset of the 1356 with reported Ea values. The prepared data has been divided into training and test sets. External test set. The test set consisted of 200 reactions randomly picked from the initial data set, to serve as a challenge for logk prediction. This included 57 structurally new transformations, that are not present in the training set. The others are reactions seen at training stage, for which the challenge consists in predicting their parameters under different conditions. Note that furthermore this test set includes 26 reactions reported at two different temperatures. The training set hence encompasses the remaining 1649 reactions. The distribution of the three properties for the training set is given on Figure 3.

Figure 3 .

 3 Figure 3. Distribution of rate constant, activation energy and pre-exponential factor values in the training set of 1649 cycloaddition reactions.

Figure 4 .

 4 Figure 4. GTM class landscape displaying separation of three different types of cycloaddition in the training set of 1649 reactions.

Figure 5 .

 5 Figure 5. GTM property landscapes for the rate constant (logk), activation energy (Ea) and the pre-exponential coefficient (logA) for the training set of 1649 reactions of cycloaddition. The encircled areas correspond to different combinations of logA and Ea contributions into logk.

  4):-(i) direct logk prediction with SVR and GTM, -(ii) prediction of logk using the Arrhenius equation (eq. 5) and the predicted values of Ea and logA,

Figure 6 .

 6 Figure 6. Predicted vs experimental logk values obtained with direct logk prediction (left) or with Arrhenius-based recalculation (right) # for the external test set of 200 reactions of cycloaddition. The statistics is given in Table3.

7 )Figure 7 .

 77 Figure 7. The comparison of the predictive accuracy for the direct GTM logk calculation and the GTM temperaturescaled landscape, for the subset of 26 reactions of cycloaddition measured at different temperatures.

Figure 8 .

 8 Figure 8. Performances of SVR and GTM models for the subset of unique structural transformations (57 reactions) obtained with the direct logk prediction method. The statistics is given in Table4.

Table 1 .

 1 Examples of reactions projected into the zones corresponding to Figure5.

Table 2 .

 2 Predictive performance of the SVR models, built separately for each of the property, and a single GTM model, universal for all three properties, in 5-fold cross validation for the training set of 1649 reactions of cycloaddition. The descriptor space is based on atom-centered fragments of the length from 2 to 3.

	Property	SVR		GTM
		R 2	RMSE	R 2	RMSE
	logk	0.94	0.45	0.78	0.86
	Ea	0.92	3.65	0.80	5.92
	logA	0.89	0.36	0.62	0.67

Table 3 .

 3 Validation of four different methods of logk assessment one test set of 200 cycloaddition reactions. 164 reactions are retained by the Fragment Control applicability domain.

	Method of logk	SVR		GTM
	assessment	R 2	RMSE	R 2	RMSE
	Direct	0.92	0.50	0.74	0.90
	Direct with AD	0.96	0.35	0.84	0.74
	Arrhenius-based	0.72	0.93	0.51	1.23
	Arrhenius-based with AD	0.93	0.51	0.83	0.79
	Temperature-scaled	-	-	0.76	0.88

Table 4 .

 4 Comparison of the direct, Arrhenius-based and temperature-scaled logk assessment for the subset of 57 unique structural transformations. Fragment Control applicability domain retained 49 reactions.

	Logk predictive	SVR		GTM
	method	R 2	RMSE	R 2	RMSE
	Direct	0.72	0.80	0.38	1.21
	Direct (AD)	0.87	0.59	0.60	1.01
	Arrhenius-based	0.33	1.26	-0.48	1.87
	Temperature-	-	-	0.40	1.15
	scaled GTM

Acknowledgement. This study was supported by Russian Science Foundation, grant No 14-43-00024. MG thanks the French Embassy in Russia for the PhD fellowship.


Carey M Lisse
email: carey.lisse@jhuapl.edu

J R Szalay
email: jszalay@princeton.edudraine@as-tro.princeton.edu

A R Poppe
email: poppe@ssl.berkeley.edu

M Horányi

V Sterken

C Beichmann

B Draine

Rosine Lallement

Anny Chantal Levasseur-Regourd

M Zemcov

M Horanyi
email: mihaly.horanyi@lasp.colorado.edu

A.-C Levesseur-Regourd

Interstellar Probe Measurements of Dust in Our Heliosphere

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Micron to sub-millimeter sized dust grains are present throughout the solar system as both interplanetary dust (IDP) and interstellar dust particles (ISD). Dust sources include grinding asteroids in the main belt, sputtered KBOs in the EKB, and sublimating comets in both regions, while sinks include solar radiation pressure and evaporation. Despite numerous previous observations over decades, the full shape and structure of the solar system's dust disks are poorly understood because we live inside of them. We especially do not understand the outer disks regions (e.g., how much dust is produced from the EKB nor how that dust migrates through the solar system) since near-Sun cometary contributions dominate near-Earth space and only one spacecraft, New Horizons, has ever flown a dedicated dust counter through the EKB. Compositional information on the makeup of interplanetary dust grains in the outer solar system is also extremely sparse, for the simple reason that when viewed from near 1 au, inner solar system grains obscure much of the signal from outer solar system dust. The ability to map the radial gradient of interplanetary dust grain composition(s) would provide strong constraints on the chemical and isotopic compositions of their parent bodies and in turn, on the formation and evolution of the solar system.

Interstellar dust grains (ISD) from the local galactic environment continuously flow through our solar system after passing through the ISM -heliosphere boundary at the heliosheath & heliopause. Their contribution to the energetics of the boundary layers are estimated to be as much as 1/3 of the total energy in the heliopause. ISDs also erode the surfaces of airless bodies at the outer reaches of our solar system, such as Oort Cloud Comets and Kuiper Belt Objects, thereby contributing to the production of the outer interplanetary dust disk. Each ISD grain also carries critical compositional information, delivering matter that may resemble the original solid building blocks of our solar system. Despite decades of observations and modeling efforts, understanding the ISD flux and its directional variability remains an unfinished and challenging task. In-situ measurements of ISD grains taken across a large range of heliocentric distances out past the Kuiper Belt are critical to understand and characterize the flux and composition of ISD, and how the solar system filters and interacts with this material. Despite this, a comprehensive experimental effort to explore the composition and dynamics of interstellar dust particles flowing through our solar system yet to be accomplished.

Using a dust analyzer instrument onboard an Interstellar Probe (ISP) to perform in situ dust collection as the spacecraft eits the solar system, we will measure the extent of the inner, near-earth zodiacal cloud; whether it connects smoothly into an outer cloud, or if there is a 2nd outer cloud sourced by the EKB and isolated by the outer planets. in situ sampling will inform about the cloud's local dust particle size and composition while performing the first ever chemical sampling of dust beyond 10 AU. It will help absolutely calibrate 3D cloud models produced by remote ISP VISIR imaging, determine if the dust in the outer system is icy or rocky or both, and help solve the current disconnect between remote imaging models of ISM dust near the Solar System and dust analyzer measurements of ISM dust made from the Ulysses and Cassini spacecraft. It will also carry a dust analyzer for the first time ever into the heliopause and heliosheath, allowing for our first understanding of the role that a dusty plasma could have in the boundary interation regions between our habitable astrosphere and the nearby local galaxy.


Carey M Lisse

M B Zemcov

A R Poppe

J R Szalay

B T Draine

M Horányi

R L Mcnutt

A Corcoros

P C Brandt

M V Paul

What can an interstellar

come

WHAT CAN AN INTERSTELLAR PROBE MISSION AT LARGE HELIOCENTRIC DISTANCES ACHIEVE WITH REMOTE IMAGING AND IN SITU DUST MEASUREMENTS?

C.M. Lisse 1 , M.B. Zemcov [START_REF] Zemcov+ | American Astronomical Society, AAS Meeting #233[END_REF] , A.R. Poppe [START_REF] Stark | [END_REF] , J.R. Szalay 4 , B.T. Draine 4 , M. Horanyi 5 , R.L. McNutt, Jr. 1 , A. Corcoros 1 , P.C. Brandt 1 , M.V. Paul 1 , V. Sterken 6 , A.C. Levasseur-Regourd 7 , R. Lallement 8 , C.A. Beichman 9 , 1 JHU-APL carey.lisse@jhuapl.edu [START_REF] Zemcov+ | American Astronomical Society, AAS Meeting #233[END_REF] RIT College of Science [START_REF] Stark | [END_REF] Space Sciences Laboratory, Univ. of California Berkeley 4 Princeton University 5 Laboratory for Atmospheric & Space Physics and Dept. of Physics, Univ. of Colorado 6 Univ. Bern 7 Sorbonne Université 8 Observatoire de Paris 9 NASA Exoplanet Science Institute, CalTech

Introduction.

Here we present the scientific potential of using an Interstellar Probe (ISP) telescope traveling to > 500 AU from the Sun to observe the brightness, shape, and composition of the dust in our solar system's debris disks and to detect, for the first time, ISM dust interacting in the heliosheath and heliopause. Discussion. Planetesimal belts and dusty debris disks are known as the "signposts of planet formation" in exosystems. The overall brightness of a disk provides information on the amount of sourcing planetesimal material, while asymmetries in the shape of the disk can be used to search for perturbing planets. The solar system is known to house two such belts, the inner Jupiter Family Comet (JFC) + Asteroid belt and the outer Edgeworth-Kuiper Belt (EKB), and at least one debris cloud, the Zodiacal Cloud, sourced by planetesimal collisions and comet evaporative sublimation.

However, these are poorly understood in toto because we live inside of them. E.g., it is not understood well how much dust is produced from the EKB since the near-Sun comet contributions dominate the inner cloud and only one s/c, New Horizons, has ever flown a dust counter through the EKB. New estimates from the NH results put the EKB disk mass at at 30 -40 times the inner disk mass [1]. Better understanding how much dust is produced in the EKB will improve our estimates of the total number of bodies in the belt, especially the smallest ones, and their dynamical collisional state. Even for the innermost Zodiacal cloud, questions remain concerning its overall shape and orientation with respect to the ecliptic and invariable planes of the solar system -they are not explainable from perturbations caused by the known planets alone.

Imaging Studies. Using new technologies and passively cooled detectors, a suitable low size, weight and power system VISNIR spectrometer/FIR imager + 10 cm class primary has been specified using a CubeSat study baseline design [2]. The VISNIR spectrometer could provide maps of the cloud's dust particle size and composition, while FIR imagery would map the dust cloud's density. 3-D cloud Mapping would occur during flythrough via tomographic inversion, and via lookback imaging once the s/c is beyond 200 AU. The lookback imaging will allow ISP to measure for the 1 st time in history the entire extent of the Zodiacal Cloud, and determine whether its inner JFC/asteroidal & outer KBO parts connect smoothly, as predicted by Stark & Kuchner [3-4] and detected by Piquette, Poppe et al. [5-8] (Figs. 23). This would also allow direct comparison of the solar system's debris disks with those observed around other nearby stars, and test theories that suggest that our solar system is planet rich but dust-poor [9]. Observing at high phase angle by looking back towards the Sun from >500 AU, we will be able to perform deep searches for the presence of rings and dust clouds around discrete sources, and thus we will be able to search for possible strong individual sources of the debris clouds -like Planet X, the Haumea family of icy collisional fragments, the rings of the Centaur Chariklo, or dust emitted from spallation off the larger KBOs. The same remote sensors will be used to map the surfaces of KBOs encountered along the way. Large-scale structure determination of the cloud should help inform us of ancient events like planetary migration and planetesimal scattering (as in the LHB), and measurement of the cloud's total brightness will allow improved removal of its signal in near-Earth cosmological measurements looking out into the Universe.

First Ever Outer Solar System

In Situ Dust Characterization. ISP can also carry the first ever in situ dust chemical analyzer past the orbit of Saturn. Based on the Europa Clipper SUDA instrument [10], it will compositionally and directionally characterize the solar system's dust clouds and will help isolate their sources, like the rocky asteroidal dust bands and the icy Haumea family fragments. Using measured dust particle masses and velocities, dust input & loss rates from these sources will be derived. Direct dust sampling will return the first ever in situ chemical analysis of EKB dust, the first ever in situ sampling of dust beyond 200 AU, and provide calibrated ground truth for cloud models produced from our imagery. It should also resolve the tension between the expected makeup of inflowing ISM microdust as determined by remote sensing and the mesasured ISM dust component found at Jupiter and Saturn by Galileo, Ulysses, and Cassini [11].

Understanding a G2V's Astropause ISM Bowshock. The role dust plays in shaping and energizing our solar system's boundary with the local galactic medium is almost completely unknown. Estimates range up to 1/3 of the energy density in the heliopause and heliosheath to be in the dust. Current models of the heliopause & sheath do not allow for the physics of a dusty plasma because the dust component is so poorly known. We do know that submicron sized dust is streaming into the solar system from the ram direction the solar system is taking through the local ISM, and the discrepancy between remote sensing models of local ISM dust and ISM dust measured inside the solar system suggests a large amount of energy is involved in diverting much of the impinging dust around the edges of the solar system in the heliosheath. Conclusions. The expected scientific return of measuring the solar system's circumsolar dust out past 500AU is large. The contributions from the Edgeworth-Kuiper Belt and the Oort Cloud, normally obscured from the Earth by locally dominant JFC and asteroid belt contributions, will be imaged and sampled. The galactic and exragalactic background VISIR fluxes will be measured clealy, without foreground contamination, for the 1 st time, as will the role of dust in the heliosheath.

Figure 1 -

 1 Figure 1 -Interstellar Probe Explorer payload at its design goal location of 1000 AU with respect to the planets, the heliopause, Alpha Centauri, and the Oort Cloud.

Figure 2 -

 2 Figure 2 -Predicted dust cloud morphologies arising from solar system JFC (JFC) & Oort Cloud (OCC) comets & Kuiper Belt (EKB) sources. (Top) Looking down on the solar system. (Bottom) Looking through the plane of the solar system. After [1].

Figure 3 -

 3 Figure 3 -In situ measurements (black data pts) and predicted dust flux contributions (colored curves) for the solar system's debris disks [1,8]. Some care should be taken in interpreting this plot, as the PIONEER 10 fluxes are for 1-100 um grains while the NH grains are 0.2 -10.0 um in size. Regardless, the overall relative shapes scale well and the predicted crossover at ~10 AU from JFC dominated to EKB dominated is seen. ISP will help us determine if another crossover from EKB dominated to OCC dominated occurs at ~100 AU, and if the EKB dust is ice, rock, & organics rich like KBOs and comets.

Figure 4 -

 4 Figure 4 -Disconnect between the nearby ISM dust size distribution predicted from remote sensing measurements (blue) and ISM dust counts measured inside the solar system (red). After [11].