Raquel Almeida

Ewa Kijak

Simon Malinowski

Zenilton K G Patrocínio

Arnaldo A Araújo

Silvio J F Guimarães

Graph-based image gradients aggregated with random forests

Gradient methods subject images to a series of operations to enhance some characteristics and facilitate image analysis, usually the contours of large objects. We argue that a gradient must show other characteristics, such as minor components and large uniform regions, particularly for the image segmentation task where subjective concepts such as region coherence and similarity are hard to interpret from the pixel information. This work extends the formalism of a previously proposed graph-based image gradient method that uses edge-weighted graphs aggregated with Random Forest (RF) to create descriptive gradients. We aim to explore more extensive input image areas and make changes driven by the RF mechanics. We evaluated the proposals on the edge and segmentation tasks, analyzing the gradient characteristics that most impacted the final segmentation. The experiments indicated that sharp thick contours are crucial, whereas fuzzy maps yielded the worst results even when created from deep methods with more precise edge maps. Also, we analyzed how uniform regions and small details impacted the final segmentation. Statistical analysis on the segmentation task demonstrated that the gradients created by the proposed are significantly better than most of the best edge maps methods and validated our original choices of attributes.

Introduction

Image segmentation is a computer vision task aiming to group pixels into regions. It is a complex task based on abstract concepts such as region coherence, perceptual similarity, and composition of objects and scenes by aggregated parts. We could define the region coherence in terms of uniformity, continuity, and contrasts that create well-defined boundaries between two adjacent regions [START_REF] Domínguez | Image segmentation: advances[END_REF]. Perceptual similarity and composition are subjective, but ideally, they would reflect our perception by hierarchically recognizing parts, objects, and concepts [START_REF] Kurzweil | How to create a mind: The secret of human thought revealed[END_REF]. This ideal could be why hierarchical methods for segmentation have remained popular since their creation [START_REF] Beucher | Watershed, hierarchical segmentation and waterfall algorithm[END_REF][START_REF] Najman | Geodesic saliency of watershed contours and hierarchical segmentation[END_REF][START_REF] Krishnammal | Wavelets and convolutional neural networks-based automatic segmentation and prediction of MRI brain images[END_REF][START_REF] Makrogiannis | A system for spatio-temporal cell detection and segmentation in time-lapse microscopy[END_REF]. Hierarchical methods build in their structure different detail levels, providing easy navigation and merging operations to build more semantically significant objects from lower-level instances. The most well-known hierarchical method is the hierarchical watershed [START_REF] Beucher | Watershed, hierarchical segmentation and waterfall algorithm[END_REF], which extends the morphological watershed [START_REF] Beucher | Use of watersheds in contour detection[END_REF] and creates a sequence of segmentations from image magnitudes.

Traditionally, the magnitudes are extracted by measuring dissimilarities in the color space (e.g., RGB, Lab) [START_REF] Perret | Evaluation of hierarchical watersheds[END_REF][START_REF] Meyer | Stochastic watershed hierarchies[END_REF] or in the gray-scale representation of the images [START_REF] Arbeláez | Contour detection and hierarchical image segmentation[END_REF][START_REF] Bosilj | Partition and inclusion hierarchies of images: A comprehensive survey[END_REF]. However, these approaches may produce a significant variation of absolute values, making it hard to determine which values compose a coherent region [START_REF] Najman | Mathematical morphology: from theory to applications[END_REF]. An alternative is to preprocess the images by applying gradient operators to emphasize desirable characteristics, such as the object's contours. Wellknown gradient operators based on kernel filters for local variation, such as Laplacian and Sobel, used to be the preferable operators [START_REF] Haris | Hybrid image segmentation using watersheds and fast region merging[END_REF][START_REF] Shanthakumar | Computer aided brain tumor detection system using watershed segmentation techniques[END_REF]. However, learned edge maps, as in the Structured Edge Detection (SED) method [START_REF] Dóllar | Fast edge detection using structured forests[END_REF], were proven to produce better results [START_REF] Perret | Evaluation of hierarchical watersheds[END_REF] and became more popular [START_REF] Perret | Removing non-significant regions in hierarchical clustering and segmentation[END_REF][START_REF] Otiniano-Rodríguez | Hierarchy based salient regions: A region detector based on hierarchies of partitions[END_REF].

Our previous work [START_REF] Almeida | Descriptive image gradient from edge-weighted image graph and random forests[END_REF] argued that although the contours are essential for contrasting adjacent regions, other characteristics reflecting uniformity and continuity are desirable for the segmentation task. We proposed a graph-based gradient operator that produced gradients with firm contours of the objects and other characteristics such as minor components, textures, and large uniform regions. We demonstrated that these gradients used as input for the watershed method produced better segmentation results than Laplace, Sobel, and SED.

Our graph-based image gradient (GIG) method uses edge-weighted graphs aggregated with Random Forest (RF) [START_REF] Breiman | Random forests[END_REF] as an image gradient operator trained on the edge detection task. The motivation behind this approach is four-fold: (i) exploit the spatial image domain preserved in the form of a grid graph: typically, we define graphs on the image domain with a structured grid adjacency relation, and the set of elemental graph components, called vertices, represents the pixels of the image; (ii) strengthen the representation with the graph relational features: we represent the adjacency relation by edges connecting neighboring vertices, and a weighting function could be associated to represent local variation. Therefore, the graph operator contains the structured spatial information and measurements of relational dissimilarities; (iii) suppress noise and encourage large consistent regions with RF: the RF algorithm, through attribute selection and implicit regularization [START_REF] Wyner | Explaining the success of AdaBoost and random forests as interpolating classifiers[END_REF], can reinforce desirable characteristics and mitigate some strong responses created by the local difference measurements on the graph operator; and (iv) maintain the analyses on the discrete space: each vertex description represents an entire neighboring region associated with a single label; therefore, we avoid training in a complex structured output space, like in [START_REF] Dóllar | Fast edge detection using structured forests[END_REF].

Our motivation goes beyond a good performance in a task. It also relies on a proposal of a machine learning framework operating on graphs, a topic of interest due to its capacity to represent multivariate information and the possibility of multiple applications, such as classification or clustering. A significant challenge is that most machine learning algorithms require regular inputs to operate, including RF. This requirement is inherently opposed to the unconstrained nature of graphs (e.g., no clear beginning or end, connected vertices are not necessarily close). The GIG representation considered the type of graph, its proximity to the original data, and the expected results, allowing us to process the graphs as regular data with a fast machine learning algorithm and avoid the long computations of graph networks.

In this work, we extend the formalism to exploit better the relationships modeled by graphs, mainly focusing on the RF mechanics and limitations. Namely, we propose:

1. Region adjacency graphs: Extends the formalism from the bijective correspondence of vertices and pixels in GIG to vertices and a set of regions produced by an initial segmentation into image super-pixels. This approach could reduce the number of data points during training and impact the gradient and the computational cost; 2. Positional features: Vertices corresponding to the pixels on the image's border have an incomplete set of neighbors. In GIG, they received padding values that disregarded the missing value's position on the regular representation. In the positional feature approach, we only take the vertices with a complete set of adjacent vertices to avoid changing the feature connotation during training; 3. Unique paths: The regular representation of the grid graph in GIG is redundant, meaning not all values are unique as the vertices on the grid path share some neighbors. Our unique path approach considers only the first instance of a neighboring vertex in a region, reducing the representation's size and allowing the region's expansion. Besides the extended formalism, we add evaluations on edge detection and comparisons with deep learning approaches [START_REF] Xie | Holistically-nested edge detection[END_REF][START_REF] Liu | Richer convolutional features for edge detection[END_REF] on edge detection and region segmentation. We aim to advocate for our assertion that a good gradient for image segmentation should present more than precise object contours by comparing it with even more accurate edge maps.

We organized this work as follows: Section 2 presents the theoretical background and Section 3 the methodology. The experimental setup and results are in Section 4, followed by a discussion in Section 5 and the conclusions in Section 6.

Theoretical background

An undirected graph G = (V, E) consists of a finite non-empty set of vertices, denoted by V , and a finite set of edges E = {{u, v} | u, v ∈ V } ⊆ V ×V . The set E induces a unique adjacency relation Γ on V , which associates u ∈ V with Γ(u) = {v ∈ V |(u, v) ∈ E}. In the image graph, the set of vertices represents the pixels of the image, and the edges represent the connections between them. Usually, we define the adjacency relation Γ by a structured component in a grid form, such as 4or 8-adjacency. We define vertex attributes as vertex functions f : V → R that map low-level image features to a vertex v.

An edge-weighted graph is denoted by (V, E, F), in which F : V × V → R is a dissimilarity function that weights the edges of G. We could assign multiple functions to the set of edges F and vertices f . This work uses F euc (u, v) = (f g (u) -f g (v)) 2 , where the vertex function f g maps the gray-scale magnitudes. Also, we apply the descriptor proposed in [START_REF] Dóllar | The fastest pedestrian detector in the west[END_REF] as a vertex function f to map the pixel's color space and color gradients.

The edge weights may represent the local variation around a vertex and serve as an image gradient operator bounded by the adjacency relation. The function defined for F euc may have a strong response to noise, as in the case of many spatial filters based on local dif- Function getAttributes(v):

1 if isTrainSet then 2 X = ones[p + 1] 3 Y = getLabel(v) 4 X[p + 1] ← Y // at
1 for vertex v in V do 2 Xv = getAttributes(v) 3 X ← append(Xv) 4 end 5 return X
ferences. The RF acting as a regularizer can diminish the noise, mitigate any eventual poor topology choice and accentuate strong connections. RF [START_REF] Breiman | Random forests[END_REF] is a non-parametric ensemble method for supervised classification and regression that relies on randomizing selected data and features and has extensive practical uses in many domains. Some authors [START_REF] Wyner | Explaining the success of AdaBoost and random forests as interpolating classifiers[END_REF][START_REF] Biau | Consistency of random forests and other averaging classifiers[END_REF] believe that randomness performs as an implicit regularization process, promoting consistency [START_REF] Biau | Consistency of random forests and other averaging classifiers[END_REF] and noise suppression [START_REF] Wyner | Explaining the success of AdaBoost and random forests as interpolating classifiers[END_REF]. The RF model consists of randomized independent trees, each trained with bootstrap samples of the input labeled data D. Applying RF to edgeweighted graph representation implies deriving regular inputs from the graph, which we will detail in Sec. 3.

Graph-based image gradients

The main challenge in this framework concerns the strategy to parse the data and create the RF's regular input without losing too much information. In our previously proposed representation, GIG [START_REF] Almeida | Descriptive image gradient from edge-weighted image graph and random forests[END_REF], the strategy depicted each vertex v ∈ V of the edge-weighted graph G = (V, E, F) as a vector X v with dimension p = |G att |, where G att is a set of selected attributes. The selection belonged to two categories: (i) vertex attributes (X f ), representing the vertices functions; and (ii) edge weights (X F ), representing the weight value in every edge v and any adjacent u ∈ Γ(v). Thus,

G att = {X f X F }.
Algorithm 1 presents the steps to create the regular GIG representation X G for an edge-weighted graph G. We repeat the procedure for all graphs in the training set and concatenate the X G outputs to make the RF training input D. Therefore, D = ((X 1 , Y 1 ), . . . , (X M , Y M )), where M is the total number of vertices in the training set. We train the RF on the edge detection task: Y ∈ {0, 1}, and all M entries have a unique discrete label. In the test step, we make the regular representation for each graph in the validation/test set and individually subject them to the estimations of the RF. The final estimated values are mapped back to the image coordinates as gradients.

We now extend the GIG formalism to propose strategies to explore more extensive input image areas without adding redundancy and make changes driven by the RF mechanics.

GIG defined on region adjacency graphs

The proposed region adjacency graph approach (GIG-RAG) reduces the number of vertices by presenting regions of grouped pixels in the set of vertices instead of a single pixel as in GIG. The GIG-RAG requires a strategy to group the pixels and considerations for edges, weights, and label attribution. Given an edge-weighted grid graph G I = (V, E, F) and a list of grouped pixels into regions S I = {r 1 , . . . , r R } for an image I, the edge-weighted RAG graph G rag has one vertex for each labeled region in S I , thus V rag = {v l | l ∈ {1, R}}. There is an edge between two vertices in V rag if an edge connects two vertices in the original grid graph G I , hence:

E rag = {{v p , v q } | v p , v q ∈ V rag ∧ p = q ∧ ∃{u, v} ∈ E | u ∈ r p ∧ v ∈ r q }. The set E rag induces a unique adja- cency relation on V rag , which associates u q ∈ V rag with Γ(u q ) = {v p ∈ V rag |(v p , u q ) ∈ E rag }.
For the weighting function, we average the edges' weights in

G I , F rag (u q , v q ) = mean{F(u, v) | ∀{u, v} ∈ E | u ∈ r p ∧ v ∈ r q }.
Finally, for label attribution in GIG-RAG, we take a majority vote of the pixels within a region to determine the region label. Also, because we have multiple vertices within each region, we do not use the vertex attributes X f ; thus, GIG-RAG attributes are only the edge weights X Frag of a regular number of closest adjacent regions.

Positional features

RF is an ensemble of multiple decision trees in which each independent tree takes local decisions to split the data considering a combination of features. The position of the features is an essential factor, as the model would assume that any subsequent data in that specific position would represent the same feature. In GIG, we added padding values to the vertices that do not have all neighbors in the grid path, disregarding the position that the missing value would assume if present. In GIG-Positional, we will take only the vertices with a complete set of adjacent vertices, creating a more regular representation for training.

Unique path

We tackle the redundancy created when we transpose the edge-weighted graph to a regular representation for our final proposal. In Algorithm 1, when we take the first and second adjacent neighbors, we inevitably cast repeated values to the regular representation as many vertices share some neighbors within the grid path. While this redundancy is not necessarily a problem, out of the 64 values obtained with two levels of neighbors with an 8-adjacency relation, only 24 of these values are unique. In the GIG-Unique approach, we will consider only the first instance of a neighboring vertex within a region. Removing the redundant values allows the expansion of the region of analysis while maintaining a similar-sized representation.

We propose three variations with an 8-adjacency relation: (i) two levels of neighbors, a 5 × 5 grid on the original image, 24 values instead of 64; (ii) three levels of neighbors, a 7 × 7 grid, 48 values instead of 512; and (iii) four levels of neighbors, a 9 × 9 grid, 80 values instead of 4096. For all variations, the regular representation has the same attributes as GIG, G att = {X V X F }, but different dimensions for the edge's attributes.

Experiments

We performed the experiments in three stages (illustrated in Fig. 1): (i) create the edge-weighted graph gradient operator from the input image; (ii) train the RF on the edge detection task to obtain the gradients; and (iii) evaluate the quality of the gradients on the segmentation task.

We chose the Berkeley Segmentation Dataset and Benchmark (BSDS500) [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF] because it proposes edge detection and segmentation labels. It contains 500 RGB images (200 train, 100 validation, and 200 test) of the same size. Each image has multiple labels performed by different annotators; thus, we applied a majority vote on the edge task to get a single label.

We evaluated GIG1 , GIG-RAG, GIG-Positional, and GIG-Unique, all with an 8-adjacency relation and without a combination of methods. The low-level descriptor discussed in Sec. 3 has 13 dimensions (|X V | = 13), therefore, for GIG and GIG-Positional p = 77. For the GIG-Unique representation, we use the terminology GIG-X -Unique, with X in {24, 80} indicating the number of unique values (thus |X F |), leading therefore to p = {37, 93}, respectively. For GIG-RAG, p = 64 (no vertex attributes), and we use the terminology GIG-RAG-R to indicate the number of desired regions of grouped pixels, where R is in {1k, 5k, 10k, 50k}. For completeness, we also evaluate GIG-Edge, the GIG variant considering only edge attributes, with p = 64.

We propose an initial segmentation into image super-pixels to group the pixels in the GIG-RAG strategy. We tested with a well-consolidated superpixel method called simple linear iterative clustering (SLIC) [START_REF] Achanta | Slic superpixels[END_REF]. SLIC is an iterative method, which clusters the pixels with the closest center, initially distributed in a regular grid, evaluates the similarity within a cluster, and recalculates the centers until convergence. SLIC execution is fast, easy to set the parameters and the number of produced regions is the closest approximation to the number of desired regions passed as a parameter.

In order to evaluate the impact of the super-pixel quality in the GIG-RAG strategy, we added some comparisons with a more modern super-pixel method: the superpixel segmentation with fully convolutional networks (SpixelFCN) [START_REF] Yang | Superpixel segmentation with fully convolutional networks[END_REF]. The SpixelFCN is a deep network trained to assign each pixel of an image, initially partitioned into a regular grid, to one of its neighboring grids. The network is an auto-encoder, in which the encoder learns the features, and the decoder aims to group pixels with similar features and enforce compactness. SpixelFCN is easy to set and overall produces better regions, but as in many deep network methods, it is limited in the number of created regions since the regular grid has a fixed size. To increase the number of super-pixels, one should increase the scale of the input image.

For the RF, we used the Random Forest Regressor included in the scikit-learn Python package [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF], which provides a parallelized implementation over the trees. For the parameters, we performed a grid search on the number of trees, the bootstrap sample size, and the number of sampled features for the split, using the validation set and the original GIG representation. The selected parameters after our search are #trees = 500, The final step in our experimental pipeline is to evaluate the quality of the gradients. We propose to apply the gradients as input to the hierarchical watershed method (HWS) [START_REF] Cousty | Incremental algorithm for hierarchical minimum spanning forests and saliency of watershed cuts[END_REF], whose performance is dependent on the gradient input (discussion in Sec. 5). Because the watershed is hierarchical, we must pass the desired number of regions to obtain the final segmentation. We used the Higra library [START_REF] Perret | Higra: Hierarchical graph analysis[END_REF] for the watershed implementation.

We present our quantitative analysis for edge detection and segmentation tasks, using the F -measure for the precision-recall metrics for boundaries and regions. We also show the Probabilistic Rand Index (PRI) [START_REF] Pont-Tuset | Measures and metameasures for the supervised evaluation of image segmentation[END_REF] measure, a metric that ponders the multiple groundtruths for the segmentation interpretation. Results are presented in terms of the optimal dataset scale (ODS), optimal image scale (OIS), and average precision (AP) through all scales. The scales for the boundaries are different thresholds applied to the edge maps to create a binary image and, for the regions, the desired number of segmented regions.

Results on edge detection

First, we want to evaluate each proposed strategy's impact on the representation performance regarding the quality of the edge maps (illustrated in Fig. 2(f)-(k)). Table 1 presents the relevant results in the validation set for the edge detection task. We omitted some similar values in Tables 1-3 to avoid repetition. All variations of the GIG-X-Unique presented similar gradients (thin contours and discreet textures) and results (less than 1.5% difference in all metrics). Moreover, all the GIG-Unique strategies had similar score metrics to the original GIG, indicating that there is not much gain in expanding the region of analysis and that the GIG redundancy does not compromise the RF generalization. The GIG-Positional results indicated that the feature position is, in fact, an essential factor during training, Table 1: F -score for boundaries in terms of optimal dataset scale (ODS), optimal image scale (OIS) and average precision (AP) through all scales (perfect scores=1). Executed on the validation set. and the edge maps created have the best performance on the task among all the compared proposals. GIG-Edge shares the number of regions with GIG while not considering the vertex attributes as GIG-RAG, and the results showed the importance of vertex attributes in the gradients and the score. We observed that the GIG-RAG strategy considerably reduced the training time proportionally to the number of regions. Nevertheless, it also reduced the performance of the edge detection task. As shown in Table 1, SpixelFCN has a slight advantage over SLIC super-pixels, but we were limited on computational resources to create a larger number of regions. For instance, to make 50k regions with SpixelFCN, one must work with images scaled to ∼ 9 times the original size. As illustrated in Fig. 3, the computed gradients vary with the number of regions, fewer regions create larger super-pixels, and the predicted label is applied to a larger area creating gradients of regions instead of contours. Increasing the number of regions creates more . GIG-RAG gradient with 50k regions computed with SLIC presents a contouroriented gradient with little texture information (i), the same is true for the GIG-Edge variation (j). GIG-Positional gradients (k) have slightly stronger borders for the main objects than GIG. In (l), the GIG gradient filtered by the morphological opening operation presenting thicker contours.

contour-oriented gradients, which is crucial for the edge detection task.

Table 2 presents the best of our proposed methods, GIG and GIG-Positional. We compared with some leading deep methods on the task: the Holistically-Nested Edge Detection [START_REF] Xie | Holistically-nested edge detection[END_REF] (HED, illustrated in Fig. 2(e)) and the Richer Convolutional Features for Edge Detection [START_REF] Liu | Richer convolutional features for edge detection[END_REF] (RCF, illustrated in Fig. 2(f)). Both HED and RCF produced better edge maps (Fmeasures reported by the authors). Still, they required considerably longer training times (measured using a GPU during 10, 000 iterations on the required augmented dataset, following [START_REF] Xie | Holistically-nested edge detection[END_REF]).

We also included results from the Structured Edge Detection method [START_REF] Dóllar | Fast edge detection using structured forests[END_REF] (SED, illustrated in Fig. 2(d)), which formalism is parallel to ours. Briefly, SED expands the RF formalism for image processing. However, to map the similarity in the structured labels of patches in the image, SED creates an intermediary reduced space, reducing the computational cost by avoiding the calculation of a continuous variance. By contrast, we made structured inputs with the image graphs to attribute a single discrete label. SED, GIG, and GIG-Positional were all trained using the same CPU, with parallelized computation over the trees in 8 CPU cores, wherein the SED's RF is composed of 8 trees and ours of 500. The training time reflects the gain of the structure input on GIG instead of the structured output on SED. Furthermore, the GIG-Positional scores were comparable to SED. For the inference time, all methods took only a fraction of a second for each image, whereas RCF and ours were slightly faster.

Results on segmentation

We evaluated the segmentation task for all the proposed GIG variants, the deep methods, and SED. In Table 3, we can see that the worst metrics are for the GIG-RAG with a small number of regions (1k and 5k) computed from SLIC. The GIG-RAG gradients computed from SpixelFCN and SLIC with a larger number of regions (10k and 50k) perform like some of the best methods on edge detection (SED and GIG-Positional). HED and GIG-Edge have similar performance on the task, whereas GIG-Edge has an advantage on the PRI metric and the AP. GIG performs better than both in all metrics. GIG-Unique under-performs compared to GIG, except for the F -measure dataset scale, meaning that we could choose a certain number of segmented regions and have more consistent results.

Finally, we have the RCF results that outperformed GIG and the others proposed in all metrics. The gradients produced by GIG and RCF are very different, from the level of details to the thickness of the contours. To investigate the thickness factor, we applied the operation opening-a well-known mathematical morphology filtering operation, consisting of one erosion to remove small regions followed by one dilation to increase object boundaries-in GIG gradients to expand the contours. We used the erosion operation with a kernel 3 × 3 to avoid enlarging small points followed by a 4×4 dilation kernel. We illustrate the result in Figure 2(l), where the GIG gradient presented thicker contours while retaining most of its details. As shown in Table 3, this operation resulted in better segmentation, indicating that thick contours are crucial to the task. Also, despite the RCF results remaining generally better, the GIG Table 2: Comparison of the best proposed methods on the edge detection task with some of the well-acknowledged methods on the dataset. We present the F -scores for boundaries in terms of optimal dataset scale (ODS), optimal image scale (OIS) and average precision (AP) through all scales, the training and inference time (per image) for all compared methods. Perfect scores=1. Executed on the test set. filtered representation outperforms RCF in the AP. It is comparable or better in the PRI metrics, which ponder areas without consent among the annotators, such as the small details better captured in GIG. We illustrate the segmentations with highlights on some critical areas in Fig. 4.

We performed statistical analysis for GIG to validate the segmentation results and present in Fig. 5 scatter graphics to illustrate. GIG is better than the best edge maps methods (SED, HED, and GIG-Positional) with statistical significance (p-values < 10e -17) and comparable to the GIG-Unique representations (p-values ∼0.02, not depicted). GIG is statistically better than GIG-Edge (p-value < 10e -14) despite both presenting similar segmentation metrics. We present the RCF comparison with the GIG filtered version, which was still inferior to the RCF (p-value < 10e -9) despite the improvement from GIG.

Discussion

Edge maps as gradients are commonly used as a preprocessing step in many applications because they are fast to compute and usually facilitate image analysis. Knowing the application and the type of analysis, one should always consider if the contour-oriented image simplification is enough for the task.

Our application is the hierarchical watershed (HWS), which deals with non-linear image analysis based on a hierarchy of partitions. In essence, given an input image (usually a gradient), it creates a multiscale image representation as an arrangement of coarse to fine partitions as an output. A ranked set of markers is necessary to construct the partitions, usually by identifying local minima based on a given geometric criterion (e.g., volume, area). We could define the hierarchy construction as the minimum spanning problem, aiming to find a solution with the minimum total cost by iteratively merging marked regions linked by the cheapest cost. The solution is optimal and consistent with the original image (details in [START_REF] Cousty | Incremental algorithm for hierarchical minimum spanning forests and saliency of watershed cuts[END_REF]).

It is important to notice that the core of the HWS resides on the cost value and the markers representing the topology of a region, both extracted from the image magnitudes. From this perspective, one could argue that the success of this method relies on a good gradient image that reflects the distribution of the orig-inal image. The usual gradients with well-delineated contours provide clear extreme values for ranking, but we argue that one should consider that these values are constantly contrasted with neighboring regions. Therefore, the depiction of uniformity and small details in conjunction with strong contours could provide additional context.

The results showed that better edge maps do not necessarily translate to better segmentation. For instance, gradients with fuzzy contours like those produced by SED and HED are not the best candidates, despite having good metrics on the edges. The small GIG-RAG representations are primarily gradients of regions instead of edges, yet, their segmentation results are not as bad as one could imagine except for the precision metric. HED and GIG-Edge have similar segmentation met-Table 3: Segmentation results for all compared methods when applied as gradient input to the hierarchical watershed method. Results presented as F -scores for boundaries in terms of optimal dataset scale (ODS), optimal image scale (OIS), and average precision (AP) through all scales, and for the Probabilistic Rand Index (PRI). Perfect F -score and PRI=1. rics but distinct gradients: HED has fuzzy, thick contours with little details, while GIG-Edge has thin, detailed contours. GIG performed better than both in all metrics: GIG shares the thick contours but not the fuzziness with HED, and it shares the details with GIG-Edge, plus additional information about patterns. To identify which characteristics command the improvement, we could examine the GIG-Unique methods, which also performed better than HED and GIG-Edge. GIG and GIG-Unique share sharp contours, details, and information patterns, indicating that these are the factors that differ and improve from HED and GIG-Edge. In contrast, GIG performs better than GIG-Unique and is varied by the thicker contours. We could also see the importance of the details and pattern information on the different results obtained with GIG-RAG from SLIC and SpixelFCN. While we could arrive at more contour-oriented gradients with a larger number of super-pixels with SLIC, the pattern information preserved with SpixelFCN gave us superior segmentation metrics and values comparable with SED and HED, even with as little as 1k regions.

Overall, gradients with thick, sharp contours perform better. But as indicated in the qualitative analysis (Fig. 4) and the precision and PRI metrics for the filtered GIG compared to RCF, additional details and information about uniform regions positively contributed to the segmentation results regarding small objects and uniformity.

Conclusions

In this work, we presented an extended formalism for the edge-weighted image gradient operator, named GIG, initially proposed in [START_REF] Almeida | Descriptive image gradient from edge-weighted image graph and random forests[END_REF]. We proposed three strategies to explore larger image areas and make changes driven by the RF mechanics to achieve a wellconsidered learning framework operating on graphs. Experiments demonstrated that reducing the number of data points by grouping the image pixels before the graph creation reduced the training time but compromised the performance on the edge detection and segmentation tasks. Also, expanding the analysis region by removing redundancy yielded similar results to the original proposal, indicating that the initial area of study already captured the necessary information, and the redundancy did not diminish the RF generalization in this particular application. Finally, the strategy that considered the position of the features regarding the RF mechanism resulted in better results in the edge detection task.

Furthermore, we extensively analyzed the gradients in the segmentation task, contrasting them with wellknown methods for the edge detection task. The results indicated that gradients with fuzzy contours yielded the worst outcomes and that uniform regions and small objects details had a substantial impact but were not as crucial as the sharp thick contours. We performed statistical analysis for GIG on the segmentation task, which is significantly better than most of the best edge maps methods. Moreover, it validated our original selection of attributes, where GIG was superior to the grouped pixels and the positional strategies and equivalent to the extended region. For future work, we will explore enriched graphs like the ones produced by hierarchical methods.

Algorithm 1 :

 1 Regular GIG and regular GIG-RAG Input : G = (V, E, F ): an edge-weighted graph, a flag isTrainSet indicating if G is a train instance and the expected representation size p. Output: X G = {(X1, Y1), . . . , (X |V | , Y |V | )}: set of regular representations of the graph vertices attributes Xv ∈ R p and its associated labels Yv in case when G is a train instance.

Figure 1 :

 1 Figure 1: The pipeline illustrated with images created at each step. We present a challenging input image (a) with transparent materials, a body of water, object occlusion, and different objects with similar colors and patterns. The graph gradient operator (b) is an image projection of the graph. The gradient in (c) is the RF predictions for the GIG method, and the segmented image (d) is the Hierarchical watershed output with 20 regions. Best viewed in color.

  GIG-24-Unique (h) GIG-80-Unique (i) GIG-RAG-50k (j) GIG-Edge (k) GIG-Positional (l) GIG filtered

Figure 2 :

 2 Figure2: Gradient computations for the input image in (a). SED, HED, and RCF present reinforced contours of the main objects but almost no details; also, most contours are fuzzy, particularly for SED and HED. GIG computes enhanced borders for large and small objects and firm representations of textures. GIG-Unique gradients present thin contours close to the ground-truth and discreet textures (g, h). GIG-RAG gradient with 50k regions computed with SLIC presents a contouroriented gradient with little texture information (i), the same is true for the GIG-Edge variation (j). GIG-Positional gradients (k) have slightly stronger borders for the main objects than GIG. In (l), the GIG gradient filtered by the morphological opening operation presenting thicker contours.

Figure 3 :

 3 Figure 3: Illustration of the gradient computations for the GIG-RAG strategy comparing the outputs of both tested superpixel algorithms. Gradients from SLIC present more apparent emphasis on the contours, while the ones from SpixelFCN preserve more texture information.

Figure 4 :Figure 5 :

 45 Figure 4: Samples of the segmentations with 10 regions produced by the RCF, GIG and GIG filtered by the morphological operation. The black areas in the ground-truth (a) indicate regions without consent among the annotators. Overall, the RCF (b) boundaries between regions are cleaner.In GIG (c) and its filtered version (d), we have more details from the object along with some background information. Also, the filtered version is more concentrated in certain details and areas of the background. In the second row, we present some region highlights to illustrate our remarks. Best viewed in color.

  Trained using GPU and F-score as reported by the authors.

		F -measure boundaries	Train time	Inference time
	Method	ODS	OIS	AP	(hh:mm:ss)	(s/image)
	SED [15]	0.712	0.724	0.750	03:53:18	0.452
	HED [21] *	0.782	0.804	0.833	11:03:42	0.215
	RCF [22] *	0.811	0.830	0.947	10:43:25	0.141
	GIG	0.635	0.661	0.648	00:09:18	0.179
	GIG-Positional	0.720	0.748	0.739	00:11:22	0.167
	*

Gradient computation code and model available at: https: //github.com/RaquelAlmeida/GIG.git

Acknowledgment

The authors thank the Conselho Nacional de Desenvolvimento Científico e Tecnológico -CNPq -(PQ 310075/2019-0 and Grant 407242/2021-0), Fundação de Amparo a Pesquisa do Estado de Minas Gerais -FAPEMIG -(Grants PPM-00006-18), Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -CAPES -(Grant STIC-AmSUD TRANSFORM 88881. 14325 8/2017-01, grant COFECUB 88887-.191730/2018-00 and funding Finance Code 001), PUC Minas and INRIA Associate Team program under the project Learning on graph-based hierarchical methods for image and multimedia data.


Ambroise De Izarra

François-Xavier Coudert
email: fx.coudert@chimieparistech.psl.eu

Alain H Fuchs

Anne Boutin
email: anne.boutin@ens.psl.eu

†pasteur

Alchemical Osmostat for Monte Carlo Simulation: Sampling Aqueous Electrolyte Solution in Open Systems

Molecular simulations involving electrolytes are usually performed at a fixed amount of salt ions in the simulation box, reproducing macroscopic concentration. Although this statement is valid in the bulk, the concentration of an electrolyte confined in a nanoporous materials such as MOFs or zeolites is greatly affected and remains a priori unknown. The nanoporous material in equilibrium with the bulk electrolyte exchange water and ions at a given chemical potential ∆µ in the semigrand canonical ensemble that must be calibrated in order to determine the concentration in the nanoporous material. In this work, we propose an algorithm based on non-equilibrium candidate Monte Carlo (NCMC) moves to ultimately perform MC simulations in contact with a saline reservoir. First, we adapt the Widom insertion technique to calibrate the chemical potential by alchemically transmuting water molecules into ions by using NCMC moves. The chemical potential defines a Monte Carlo osmostat in the semigrand constant volume and temperature ensemble (∆µ, N, V, T) to be added in a Monte Carlo simulation where the number of ions fluctuates. In order to validate the method, we adapted the NCMC move to determine the free energy of water solvation and subsequently explore thermodynamics of electrolyte solvation at infinite dilution in water.

Finally, we implemented the osmostat in MC simulations initialized with bulk water that are driven towards electrolytes of similar concencentration as the saline reservoir.

Our results demonstrate that alchemical osmostat MC simulation is a promising tool for use to sample electrolyte insertion in nanoporous materials.

Introduction

In computational chemistry, atomistic Monte Carlo (MC) simulations are routinely used to study the structure and equilibrium properties of condensed phases, solid or liquid, and their interfaces. Water is an important fluid in nature, and aqueous electrolytes at various concentrations are a key component of many chemical and biological systems, as well as used in applications such as batteries, desalination, and pollution remediation. Consequently, the insertion of ionic species in a dense fluid such as water is required for molecular simulation in many situations involving open systems, such as the osmotic equilibrium through membranes or the uptake of electrolytes in nanoporous materials. The Grand Canonical Monte Carlo (GCMC) simulation method is well-suited for simulating the molecular insertion or exchange of neutral species in fluid phases, and is used ubiquituously in the simulation of gas and liquid adsorption in porous materials. [START_REF] Ren | Thermodynamic exploration of xenon/krypton separation based on a high-throughput screening[END_REF] However, its use for ionic species is particularly inefficient: in addition to the fact that finding suitable cavities large enough in the fluid to accomodate ion insertion would appear infrequently, as for any large molecule, the strong organization of the solvent induced by the electric charge gives the insertion move a low probability. Secondly, electrical neutrality of the system imposes the insertion of pairs of charge-balancing ions (or more), which further decreases the probability of a successful Monte Carl insertion move.

Yet, the convergence and accuracy of a MC-based molecular simulation depends greatly on the ability to carry out insertion or deletion of ions, in order to describe phenomenan like the adsorption of electrolytes at interfaces or in porous materials, without an explicit description of the electrolyte reservoir (which is computational expensive, and induces finitesize effects). The design of adapted MC moves that have higher probabilities of particle insertion or deletion in dense fluid is therefore required, 2 for which a chemical potential is imposed and dictate composition evolution in response to this chemical potential. A large number of biasing techniques have been developed to this end, and without attempting an exhaustive review of this vast literature, we summarize below some of the most relevant works to provide the reader with context for our own work. able energetic location, [START_REF] Shelley | A configuration bias Monte Carlo method for ionic solutions[END_REF][START_REF] Snurr | Prediction of adsorption of aromatic hydrocarbons in silicalite from grand canonical Monte Carlo simulations with biased insertions[END_REF] suitable cavities [START_REF] Mezei | A cavity-biased (T, V, µ) Monte Carlo method for the computer simulation of fluids[END_REF] or distance-biased insertion holes 3 in the system to enhance sampling of ion insertion in water. An important example is the method of thermodynamic integration originally proposed by Kirkwood [START_REF] Kirkwood | Statistical mechanics of fluid mixtures[END_REF] which is suited for calculations of free energies between two given states. In past work, it has been applied to the calculation of the chemical potential of ions in aqueous solutions. [START_REF] Hummer | Free energy of ionic hydration[END_REF][START_REF] Ferrario | Solubility of KF in water by molecular dynamics using the Kirkwood integration method[END_REF][START_REF] Lynden-Bell | Using simulation to study solvation in water[END_REF] Its basic principle is the definition of a thermodynamic path where an ion is gradually added to the electrolyte solution by coupling its interaction with a scaling parameter λ, varying in the interval [0,1], during a MC or molecular dynamics (MD) simulation. Other advanced techniques such as multicanonical methods, [START_REF] Berg | Multicanonical ensemble: A new approach to simulate firstorder phase transitions[END_REF][START_REF] Orkoulas | Phase behavior of the restricted primitive model and square-well fluids from Monte Carlo simulations in the grand canonical ensemble[END_REF][START_REF] Fenwick | Expanded ensemble and replica exchange methods for simulation of protein-like systems[END_REF] metadynamics for chemical potential calculation [START_REF] Perego | Chemical potential calculations in dense liquids using metadynamics[END_REF] or advanced configurational bias moves [START_REF] Siepmann | Configurational bias Monte Carlo: a new sampling scheme for flexible chains[END_REF][START_REF] Laso | Simulation of phase equilibria for chain molecules[END_REF] have been developed.

Another important class of methods are the fractional component Monte Carlo moves [START_REF] Nezbeda | Recent progress in molecular simulation of aqueous electrolytes: Force fields, chemical potentials and solubility[END_REF][START_REF] Shi | Continuous fractional component Monte Carlo: an adaptive biasing method for open system atomistic simulations[END_REF] where gradual insertion and deletion of particles are performed through use of a coupling parameter λ with an adaptative bias potential -bearing some similarity in essence to the thermodynamic integration methods. These methods have been used in single and multicomponent systems in dense phases, as a way to improve phase equilibria simulations and to allow the direct calculation of chemical potentials and partial molar properties. In fractional component MC simulations, the number of particles is a continuous variable through λ in a way that particle insertion or deletion is done with a so-called fractional particle. More precisely, a set of Metropolis-like acceptance rules can be derived for λ: if a change of the coupling parameter occurs so that λ > 1, the current fractional component particle is converted into a real particle and a new fractional particle in inserted with coupling parameter of λ-1. Likewise, if the coupling parameter becomes λ < 0, the fractional particle is removed from the system and a remaining particle is designed as the new fractional particle with a coupling parameter of λ + 1. To further improve the convergence of these simulations, an appealing strategy is to improve the acceptance probability of particle insertions or deletions by using some local relaxation process, in order to minimize energy penalty and remove non-physical forces by collective motion in the surrounding molecules. This relaxation is typically performed through the introduction of a number of MD steps. Therein, C ¸agin and Pettitt developed an hybrid method [START_REF] C ¸agin | Molecular dynamics with a variable number of molecules[END_REF] where an MC move is used to change the composition of the system, associated with deterministic MD trajectories for relaxation. Another example of scheme has been proposed by Boinepalli and Attard 21 where a MD step in the microcanonical ensemble is performed at a given state associated to a certain value of coupling parameter λ of the current fractional molecule. For further details, we refer the interested reader to the recent and comprehensive review by Rabhari and co-workers on fractional component Monte carlo methods. [START_REF] Rahbari | Recent advances in the continuous fractional component Monte Carlo methodology[END_REF] Recently, Nilmeier and co-workers [START_REF] Nilmeier | Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation[END_REF] introduced a new type of hybrid MC move denoted nonequilibrium candidate Monte Carlo (NCMC) based on non equilibrium molecular dynamics. In this method, a candidate move is proposed through a finite-time process where the system is driven out of equilibrium, and accepted with a probability that preserves microscopic reversibility. In this process, the acceptance rule follows a Metropolis-type probability, where the total work performed during the out-of-equilibrium process is taken into account (instead of the difference of free energy). The generation of out-of-equilibrium states induces an additional computational cost, compensated by the reduced structural correlation that can occur in dense systems such as electrolytes and the versatile nature of the algorithm.

Along this line, Ross and co-workers 24 adapted the NCMC scheme through implementation of a Monte Carlo osmostat added to MC simulations that can sample from a semigrandcanonical isobaric isothermal ensemble in which the number of NaCl salt pairs (that we will called later an electrolyte unit for generalization) varies to study ion distribution in water around biomolecules. The NCMC moves enable high acceptance rates for which water molecules are alchemically transformed into a pair of salt ion (i.e. an electrolyte unit) or vice-versa.

MC simulation in the semigrand canonical ensemble requires the calibration of the chemical potential ∆µ used as a thermodynamic parameter governing the osmostat in equilibrium with an infinitely sized reservoir of electrolyte at a specified salt concentration. We invite the readers to refer to the original publication of Ross and Chodera for more details. In their work, the self-adjusted mixture sampling [START_REF] Tan | Optimally adjusted mixture sampling and locally weighted histogram analysis[END_REF] sample uniformly salt concentration during a single simulation, follow by application of the Bennett acceptance ratio to extract relevant free energy difference from NCMC move. The chemical potential can be thus extracted from relation between macroscopic electrolyte concentration and free energies.

In this work, we proposed a variant of the osmostat implementation that we implemented in RASPA, [START_REF] Dubbeldam | RASPA: molecular simulation software for adsorption and diffusion in flexible nanoporous materials[END_REF] a software dedicated to simulation of adsorption in porous materials. Our ultimate goal would be to employ in the future our osmostat version for the study of various electrolyte insertion in porous materials in equilibrium with a saline reservoir. The paper is organized as follows. First, we proposed another version of chemical potential calibration in the semigrand canonical ensemble (∆µ, N, V, T ) by using the Widom insertion technique with NCMC move to transform alchemically water into an electrolyte unit. Then, in order to assess the accuracy of our chemical potential calibration, NCMC move is used to study thermodynamics of electrolyte solvation at infinite dilution by calculating the free energy of water solvation with our modified Widom method. Finally, the chemical potential is used as the thermodynamic parameter governing the osmostat to gradually perturb bulk water that transition into an electrolyte that reaches concentration dictated by the previously calibrated chemical potential.

General simulation methodology

We ran molecular dynamics (MD) and Monte Carlo (MC) simulations using the RASPA software. [START_REF] Dubbeldam | RASPA: molecular simulation software for adsorption and diffusion in flexible nanoporous materials[END_REF] We used a combination of Lennard-Jones and electrostatic potentials. Water molecules are represented with the rigid TIP4P model [START_REF] Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF] and the electrolytes are modeled as charged particles that differs by their Lennard-Jones parameters and charges. We considered the following monovalent electrolytes which are: LiCl, NaCl, KCl, RbCl, CsCl, KF, KBr and KI. We also explored the following divalent electrolytes to probe influence of charge ion:

We used Lorentz-Berthelot mixing rules for cross-terms in Lennard-Jones potential, and Ewald summation was used to account efficiently for electrostatic interactions. We used a cutoff of 8.0 Å for both the Lennard-Jones potential and the separation between real space and Fourier space in the Ewald summation. MD simulations were performed with a timestep of 1 fs and all simulations used a Nosé-Hoover thermostat, setting the temperature at 300 K.

In order to explore the phase space by using MC simulations, translations were performed both on ions and water while rotations were applied only to water molecules, in a ratio of 1:1:1 respectively. NCMC moves are performed less frequently than these "standard" MC moves so that the relative probability to perform NCMC moves was set to 10 -3 . The reason is that ion insertion/deletion achieved during a NCMC move would appear less frequently than simple MC moves to rotate or translate a chosen particle. Also, short MD steps performed during alchemical transformation of water molecules into ions (or vice versa) accounts also for collective motion of particles to accommodate insertion of water (or ions). Consequently, the execution time of a NCMC move is much higher than the execution time of a "standard" MC move, as discussed below. In our validation of the method, we considered systems containing at most 150 molecules to keep the computational cost reasonable, due to the sequential nature of the RASPA code.

Calibration of the chemical potential in the semigrand canonical ensemble

Let us consider a system made of N molecules with N H 2 O water molecules, N C m+ salt cations and N A n-salt anions where m and n design respectively the absolute charge of the cation and of the anion. The number of electrolyte units (i.e. the number of salt pair of NaCl where m and n are equal to one for instance) is determined as follows:

N CnAm = N C m+ + N A n- m + n (1)
The partition function for such a system having N CnAm electrolyte units reads:

Z(N CnAm , N, V, T ) = N ! N C m+ !N A n-!N H 2 O ! dr e -βU [r,N CnAm ] (2)
The chemical potential in the semigrand canonical ensemble ∆µ (n+m)H 2 O→CnAm (denoted ∆µ CnAm for clarity) is the difference of chemical potential between withdrawing (n+m) water molecules and inserting an electrolyte unit C n A m . [START_REF] Ross | Biomolecular simulations under realistic macroscopic salt conditions[END_REF] The number of water molecules and electrolyte units become respectively

N ′ H 2 O = N H 2 O -(n+m) and N ′ CnAm = N CnAm +1 which
implies that the number of anions and cations increases respectively to

N ′ A n-= N A n-+ m and N ′ C m+ = N C m+ + n.
Finally, the chemical potential ∆µ CnAm is expressed as follows:

∆µ CnAm = - 1 β ln Z(N ′ CnAm , N, V, T ) Z(N CnAm , N, V, T ) ∆µ CnAm = - 1 β ln     m+n-1 i=0 (N H 2 O -i) n j=1 (N C m+ + j) m k=1 (N A n-+ k) e -βW (n+m)H 2 O→CnAm N CnAm     (3)
where e -βW (n+m)H 2 O→CnAm N CnAm is the Boltzmann factor of the work performed for transforming alchemically water molecules into an electrolyte unit (denoted alchemical work from now on), averaging in the canonical ensemble of the system containing N CnAm electrolyte units. This adapted version of the Widom technique is slightly different from the original one applied to GCMC simulations as it requires only the difference of potential energy before and after insertion of a particle. [START_REF] Widom | Some topics in the theory of fluids[END_REF] In our case, we employ NCMC moves [START_REF] Nilmeier | Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation[END_REF][START_REF] Ross | Biomolecular simulations under realistic macroscopic salt conditions[END_REF] to compute the alchemical work W (n+m)H 2 O→CnAm (denoted W CnAm for clarity) performed during the alchemical transformation: (n + m) water molecules were chosen randomly in the box and alchemically transformed into an electrolyte unit C n A m over an alchemical path divided in T = 200 segments, each one consisting of a NCMC step where the system is perturbated and then relaxed. For a given segment t, the alchemical path is decribed by two parameters: the state of the system r t and the values of the interpolated non-bonded parameters λ t of the molecules undergoing alchemical transformation. In the same manner as in the work of Ross and co-workers, [START_REF] Ross | Biomolecular simulations under realistic macroscopic salt conditions[END_REF] the alchemical path is a linear interpolation of the nonbonded parameters of water and salt ions, applied to randomly chosen (n+m) water molecules to be transformed into ions:

λ t = (1 -f t )λ water + f t λ ion and f t = t T (4)
A segment t within the alchemical path consists in a perturbation of the system by updating non-bonded parameters (r t , λ t )→(r t , λ t+1 ), followed its relaxation from (r t , λ t+1 )→(r t+1 ,

W CnAm = T t=1 U [r t , λ t+1 ] -U [r t , λ t ] (5)
Hence, a single NCMC move yields to one value of W CnAm and we performed 3000 NCMC moves to estimate W CnAm and obtain the chemical potential ∆µ CnAm during a MC simulationn by using our modified Widom technique. We present typical execution time of a NCMC move compared to standard moves implemented in RASPA in the Supplementary Information (Table S8). A NCMC move is computationally expensive as its execution time is seven order of magnitude higher than those for translation and rotational moves, at the benefit of a proper equilibration of the system during the alchemical transformation. The calibration process of the chemical potential is performed for different electrolyte over a wide range of concentration (Table 1). The Figure 1 shows the evolution of the alchemical work and chemical potential for various electrolytes as a function of the electrolyte concentration. We could not explore KI electrolyte for concentration greater than 10 M (Figure 1-B, dark blue) because the largest anions I -(σ I -= 5.167 Å) could not be inserted in the box with other species. For all electrolyte, the alchemical work and the free energy are the highest (i.e., more negative) after insertion of a single electrolyte unit (c = 0.37 M), which is consistent with similar calculations performed on NaCl electrolyte. [START_REF] Ross | Biomolecular simulations under realistic macroscopic salt conditions[END_REF] Upon increase of electrolyte units initially present in solution, the alchemical work required to transform water into ions decreases in absolute value as well as the chemical potential. Previous simulations of potassium halide aqueous solutions [START_REF] Soper | Ion solvation and water structure in potassium halide aqueous solutions[END_REF] showed that the hydration number (i.e. the number of water molecules in the first hydration shell) of cations and anions drops while the number of ion pair increases slightly. This trend is consistent with the radial distribution functions (Figures S4 to S15

in the Supporting Information) and the number of water molecules in the first solvation shell around ions (Tables S9 to S20). The hydration shells around ions are disturbed when concentration of electrolyte increases, which lower solvation energy of salt ions. Another common trend highlighted is that the alchemical work and chemical potential increases in absolute value when the radius of the ion diminishes, because smaller ions have a higher solvation free energy. [START_REF] Irudayam | Prediction and interpretation of the hydration entropies of monovalent cations and anions[END_REF] The absolute difference of chemical potential after insertion of an electrolyte unit is about ∼ 215 kJ/mol between LiCl and CsCl, and about ∼ 494 kJ/mol between MgCl 2 and BaCl 2 , which is more than twice the value for monovalent electrolytes.

We recall that Li + and Mg 2+ have similar sizes of respectively σ Li + = 1.505 Å and σ Mg 2+ = 1.630 Å, as well as Cs + and Ba 2+ with a respective size of σ Cs + = 3.883 Å and σ Ba 2+ = 3.820 Å. A higher charge of the cation leads thus to a stronger energy of solvation. This observation is in qualitative agreement with previous experiments [START_REF] Näslund | X-ray absorption spectroscopy study of the hydrogen bond network in the bulk water of aqueous solutions[END_REF][START_REF] Omta | Negligible effect of ions on the hydrogen-bond structure in liquid water[END_REF][START_REF] Collins | Ions in water: Characterizing the forces that control chemical processes and biological structure[END_REF] and simulations [START_REF] Irudayam | Prediction and interpretation of the hydration entropies of monovalent cations and anions[END_REF] showing that only the first solvation shell is perturbed by monovalent cation, which is not the case for divalent cations where perturbations occur in the second solvation shell, [START_REF] Irudayam | Prediction and interpretation of the hydration entropies of monovalent cations and anions[END_REF] that could account for the higher difference of chemical potential for divalent cation electrolytes. The evolution of the alchemical work for all alkali-chloride electrolyte solutions as a function of the NCMC moves is given in Figure 2 and in the SI for every electrolytes (Figures

S1 to S3

). The numerical data and an estimation of the error for both the alchemical work and the chemical potential are also provided in SI (Tables S4 to S7). The evolution of the alchemical work as a function of NCMC moves remains constant but displays more fluctuations when the number of electrolyte units increases in the system, leading to higher error calculated for these cases. Overall, errors of the alchemical work calculated for monovalent electrolytes are below of ∼ 1 kJ/mol, except for divalent cations where error reach at most ∼ 4 kJ/mol. Along this line, errors associated with the calculation of the chemical potential are reasonable and reach at most ∼ 5.5 kJ/mol (∼ 1.3 kcal/mol) at high electrolyte concentration. However, the error associated with the calibration process using our method is higher than those of Ross and Chodera for which error of ∆µ NaCl when N NaCl is four times lower than ours. However, our calibration process is rather straightfoward to implement and errors are small enough so that it represents another way for calculating chemical potential in solution with low electrolyte concentration.

Thermodynamics of electrolyte solvation at infinite dilution

We consider one of the previous system for which we calculated the chemical potential ∆µ CnAm where only one electrolyte unit (N CnAm = 1) is inserted in the box (first line of upper and lower Table 1). Electrolyte solvation at infinite dilution is thus studied using a moderate box size. These small systems contain enough water molecules to capture the thermodynamics of electrolyte solvation because experiments [START_REF] Näslund | X-ray absorption spectroscopy study of the hydrogen bond network in the bulk water of aqueous solutions[END_REF][START_REF] Omta | Negligible effect of ions on the hydrogen-bond structure in liquid water[END_REF][START_REF] Collins | Ions in water: Characterizing the forces that control chemical processes and biological structure[END_REF] and MD simulations [START_REF] Irudayam | Prediction and interpretation of the hydration entropies of monovalent cations and anions[END_REF] demonstrated that ions disturb hydrogen bonds between water molecules only in the first and second solvation shell. At the end of each NCMC move, we capture the distance separating the ion pair in order to construct the radial distribution function and deduce the potential of mean force of ion pair at infinite dilution, as shown in Figure S16 in the Supporting Information. This approach is more suitable than analyzing (N, V, T ) simulations where only one ion pair is present because ion pair formation (i.e., short distance between ions) can be more easily sampled.

We first calculated the free energy of electrolyte solvation in water at infinite dilution.

This quantity is the variation of the Helmholtz free energy ∆F vacuum→CnAm (abbreviated later by ∆F CnAm for clarity) for bringing an electrolyte unit C n A m from the vacuum to the TIP4P liquid water. We divided the solvation process in two steps. First, water molecules are removed from vaccum and solvated in bulk water. Then, water molecules are alchemically transformed into ions as described in the previous section, so that the free energy of solvation reads:

∆F CnAm = ∆µ CnAm + ∆µ (n+m)H 2 O -(n + m)k B T (6)
where ∆µ (n+m)H

where e -βW (n+m)H 2 O (150-(n+m))H 2 O is the Boltzmann factor of the water solvation work averaging in the canonical ensemble of the system containing 150 -(n + m) water molecules. However, the NCMC move is slightly modified as follows: the calculation of the solvation work W (n+m)H 2 O was performed in two stages, [START_REF] Grossfield | Ion solvation thermodynamics from simulation with a polarizable force field[END_REF] growth and charging. The molecules are grown according to a path of T = 200 segments where each segment t consists in a perturbation of the system by updating the Lennard-Jones parameters followed by relaxation through 200 NVE MD steps:

ϵ ij t = f t ϵ ij final σ ij t = f t σ ij final f t = t T (8)
Secondly, another path of T = 200 segments is repeated where the electrostatic interaction is turned on, in the same way as for the growth of the molecules. The water solvation work is summed during the process of growth and charging in the same way as equation 5. (or -6.334 kcal/mol), which is consistent with previous studies [START_REF] Ross | Biomolecular simulations under realistic macroscopic salt conditions[END_REF][START_REF] Jorgensen | Free energy of TIP4P water and the free energies of hydration of CH 4 and Clfrom statistical perturbation theory[END_REF] and confirms that the NCMC method is versatile as it can be used to performed both alchemical transformation and particle insertion.

The evolution of W 2H 2 O and its decomposition in the work performed during growth and charging as a function of NCMC moves is presented in Figure 3. It is clear that the variation of the non-bonded parameters between two alchemical segment is small enough to prevent blow-up of the system during the relaxation process by NVE MD as shown by the fluctuating but constant value of the work performed during growth (19.542 ± 0.184 kJ/mol) and charging (-68.309 ± 0.112 kJ/mol). Finally, one can deduced the free energy of solvation for electrolytes at infinite dilution, summarized in left column of Table 2. Overall, our calculations are in qualitative agreement with experimental results of Schmidt et al 41 in that the free energy are all negatives, and become less negative for larger ions. [START_REF] Irudayam | Prediction and interpretation of the hydration entropies of monovalent cations and anions[END_REF] Also, Ross and Chodera calculated solvation free energy of NaCl at infinite dilution in TIP4P water and found -713.79 kJ/mol [START_REF] Ross | Biomolecular simulations under realistic macroscopic salt conditions[END_REF] which are within 5% of our estimates (-733.4 KJ/mol). In addition, we estimated separately the variation of internal energy ∆U CnAm and deduced the variation of entropy ∆S CnAm from the variation of free solvation energy ∆F CnAm . The variation of free energy is the sum of the variation of internal energy and entropy

∆F CnAm = ∆U CnAm -T ∆S CnAm (9)
By identifying equation 9 with equation 6, the expression for the variation of internal energy writes 35

∆U CnAm = U CnAm+(150-(n+m)H 2 O) -U (150-(n+m)H 2 O) - 1 2 (n + m)k B T (10)
We calculated separately the internal energies of the system containing only 150-(n+m)

water molecules and with an electrolyte unit C n A m in addition to 150 -(n + m) water molecules by performing MD simulations in the NVT ensemble for 20 millions steps to gather statistics after equilibration during 1 millions steps. By taking the difference between the variation of Helmholtz free energy and variation of internal energy, we obtain an estimation of the variation of entropy ∆S CnAm for electrolyte solvation at infinite dilution. The results for the variation of internal energy and of entropy are displayed respectively in the middle and left columns of Table 2. Again, the values for ∆U CnAm and ∆S CnAm are all negative and are more negative for smaller ions, except for CsCl that displays a more negative entropy than that their hydration number differs greatly. [START_REF] Duboué-Dijon | Hydration and ion pairing in aqueous Mg2+ and Zn2+ solutions: force-field description aided by neutron scattering experiments and ab initio molecular dynamics simulations[END_REF] This is because divalent cations present strong polarization and charge transfer effects that can not be modeled with simple empirical force fields. Better force fields has been developed to model properly divalent cations in solution.

Li and Merz proposed to add another term scaling as r -4 to the Lennard-Jones potential [START_REF] Li | Parameterization of highly charged metal ions using the 12-6-4 LJ-type nonbonded model in explicit water[END_REF] and obtained a better estimation of the hydration free energy, ion oxygen distance as well as coordination number in system involving multivalent cations than using empirical force field with regular 6-12 Lennard-Jones potential. Polarizable force field were successfully applied to divalent cations for which the charge and the van der Waals centers are distributed to multiple site represented by dummy atoms. [START_REF] Saxena | Multisite ion models that improve coordination and free energy calculations in molecular dynamics simulations[END_REF] Also, the Electronic Continuum Correction is another approach that takes into account of the polarization of the divalent cation in a mean field approach: the charge of the cation are rescaled by the inverse square root of the electronic part of solvent dielectric constant. [START_REF] Leontyev | Accounting for electronic polarization in nonpolarizable force fields[END_REF] In this work, we have made the choice of using the refined empirical force field proposed by Mamatkulov and coworkers, [START_REF] Mamatkulov | Force fields for divalent cations based on singleion and ion-pair properties[END_REF] not only because it reproduces empirically experimental energy solvation at infinite dilution for divalent cation electrolytes, but also for the transferability of the force field for monovalent salt, water and later, porous materials.

MC simulation in the semigrand canonical ensemble: from bulk water to electrolyte solution

We implemented the osmostat in RASPA so that electrolyte units can be added or removed during a simulation, in the same way as in the work of Ross and co-workers. [START_REF] Ross | Biomolecular simulations under realistic macroscopic salt conditions[END_REF] By setting up the osmostat, a system is driven to chemical equilibrium with a saline reservoir at a given electrolyte concentration that required the calibration of ∆µ CnAm previously done in the first section. The maximal number of electrolyte units N max CnAm that can be present in the system is the floor function of N H 2 O /(n + m) where N H 2 O is the number of water molecules in the bulk system chosen as an initial configuration (150 water molecules). We briefly present a description of the osmostat implementation in RASPA.

• First, we choose whether to insert or delete an electrolyte unit

P insert =            1 if N CnAm = 0 1 2 if 0 < N CnAm < N max CnAm 0 if N CnAm = N max CnAm ( 11
)
of the corresponding probability densities (equation 13):

α(N CnAm → N CnAm + 1) = π(r, N CnAm + 1, N, ∆µ CnAm , V, T ) π(r, N CnAm , N, ∆µ CnAm , V, T ) α(N CnAm → N CnAm + 1) = m+n-1 i=0 (N H 2 O -i) n j=1 (N C m+ + j) m k=1 (N A n-+ k) e -β(W CnAm -∆µ CnAm ) (14)
In the same manner, the probability of acceptance of a trial NCMC move in which an electrolyte unit is transmuted into (n + m) water molecules denoted α(N CnAm →

N CnAm -1) writes α(N CnAm → N CnAm -1) = π(r, N CnAm -1, N, ∆µ CnAm , V, T ) π(r, N CnAm , N, ∆µ CnAm , V, T ) α(N CnAm → N CnAm -1) = n-1 i=0 (N C m+ -i) m-1 j=0 (N A n--j) m+n k=1 (N H 2 O + k) e -β(W -CnAm +∆µ CnAm ) (15)
Finally, the equations 14 and 15 lead to the MC moves implemented for insertion or deletion of an electrolyte unit that define the osmostat

acc(N CnAm → N CnAm + 1) = min     1, m+n-1 i=0 (N H 2 O -i) n j=1 (N C m+ + j) m k=1 (N A n-+ k) e -β(W CnAm -∆µ CnAm )     acc(N CnAm → N CnAm -1) = min      1, n-1 i=0 (N C m+ -i) m-1 j=0 (N A n--j) m+n k=1 (N H 2 O + k) e -β(W -CnAm +∆µ CnAm )      (16)
We S4 and Figure S1) and more rejection of NCMC moves. Nonetheless, our osmostat implementation in RASPA is able to efficiently sample the electrolyte concentration within a small amount of MC cycles.

Conclusion

A technique for insertion of electrolytes in water has been implemented in this work, based on non-equilibrium candidate Monte Carlo moves as proposed in the previous work of Ross and co-workers. At constant number of particles in the system, water molecules and electrolytes are exchanged at a given chemical potential, so that MC steps that insert or delete ions define an osmostat in contact with the system.

First, the calibration of the osmostat has been handled using the Widom insertion technique requiring alchemical work extracted from proposed NCMC moves where water molecules are alchemically transformed to an electrolyte unit. Our method is straightforward to implement, at the expense of a slightly higher error on the computed alchemical potential compared to the algorithm of Ross and co-workers. Then, we calculated solvation free energy of water by decoupling Lennard-Jones and electrostatic interactions during insertion of water molecules in the bulk using NCMC moves. Previous chemical potential for inserting a single electrolyte unit in bulk water by alchemical transformation is used to grasp thermodynamics of electrolyte solvation at infinite dilution. Our value of water solvation free energy as well as electrolyte solvation free energies are in good agreement with the literature, demonstrating that NCMC move is a versatile tool that can be used for both alchemical transormation and particle insertion in dense fluid. Finally, the MC moves for inserting/removing an electrolyte unit are similar to those of Ross and co-workers and show that our implementation of the osmostat yields the correct electrolyte concentration in system initialized with pure bulk water.

We chose to develop and test our own version of the osmostat in RASPA to have access of functionalities already implemented in the software for studying adsorption of electrolytes solutions in nanoporous materials. We plan to use our osmostat to simulate adsorption of various electrolyte solutions in zeolithes or MOFs that requires calibration of osmostat of different electrolytes solution over a wide range of concentration, which was the main goal of the current work. This opens the route for probing thermodynamics of electrolyte intrusion, pression intrusion and structure of intruded electrolyte solution in nanoporous material which will be the object of the next study. However, because of the choice of our empirical force field that is not optimal for divalent cations, we do not focus on reproducing precisely the structure of divalent cation electrolytes although we expect to grasp the thermodynamics and thus concentration of divalent cation electrolyte in porous materials. In addition, a version of RASPA with implementation of the osmostat will be provided therein. As a last remark, the osmostat was tested on relatively small boxes of monovalent or divalent-cation aqueous electrolytes as a proof-of-concept because RASPA is a sequential program and it would be time consuming to consider larger systems. In this regards, zeolites/MOFs are in most cases periodic so that only a supercell of several repeating units is required for future MC simulations where only few hundreds of free particles are adsorbed which is fortunately feasible with RASPA.

λ

  t+1 ) by using 200 steps of MD in the NVE ensemble. A NCMC move requires thus 40 000 MD steps which represents a duration of 40 ps. At each segment, the change of potential energy after perturbation of the system is recorded and summed, yielding to the alchemical work W CnAm performed on the system to alchemically transform (n + m) water molecules into an electrolyte unit C n A m along the alchemical path:

Figure 1 :

 1 Figure 1: A-B) For monovalent electrolytes CA: Evolution of the alchemical work performed during alchemical transformation W CA (left) and associate chemical potential ∆µ CA (right).For A), the cation is varied and the anion is fixed to Cl -. For B), the cation is fixed at K + and the anion is varied. C) For divalent electrolytes CA 2 : Evolution of the alchemical work W CA 2 (left) and associate chemical potential ∆µ CA 2 (right).

Figure 2 :

 2 Figure 2: Evolution of the alchemical work as a function of the NCMC moves for alkalichlorides electrolyte solution at two different concentration. Left) N CA = 1 is equivalent to 0.37 M. Right) N CA = 30 is equivalent to 18.45 M.

Figure 3 :

 3 Figure 3: Evolution of alchemical work W 2H 2 O (in black) for solvating 2 water molecules in bulk water. The contribution of vdw and electrostatic interactions are displayed in blue and red, respectively.

  RbCl, likely due to accumulated uncertainty from ∆U CnAm and ∆F CnAm . Not only ∆S CnAm are underestimated for all electrolytes (less negative) compared to experimental data, but also values of entropy for divalent cation electrolytes show a larger difference with experiments than monovalent electrolytes which is likely due to the limitation of the empirical force field used in this work. It is well known that empirical force fields that model ions with their charge and a Lennard-Jones center do not describe properly divalent cations in solution. For instance, Mg 2+ and Zn 2+ have almost the same radius (0.72 vs 0.74 Å), but it has been shown

Figure 4 :

 4 Figure 4: Evolution of CsCl concentration for different chemical potential ∆µ CsCl as a function of MC cycles, for which the concentration of the saline reservoir is indicated in black for each case. The distribution of salt equilibrium is given on the right of the time-series plots.

Table 1 :

 1 Composition of systems of 150 particules after transformation of water into electrolyte units. For each concentration, ∆µ CnAm is calculated according to equation 3. The upper table presents concentration of monovalent electrolytes (n = 1 and m = 1) while the lower table presents concentration of divalent cation electrolytes (n = 1 and m = 2).Number of electrolyte units (NCA ) N H 2 O ratio N H 2 O /N CA Concentration (M) CA 2 ) N H 2 O ratio N H 2 O /N CA 2 Concentration (M)CnAm and the number of water N H 2 O shown in Table1are the targeted number of species after removal of (n + m) water and insertion of an electrolyte unit C n A m .

	1	148	148	0.37
	6	138	23	2.41
	11	128	11.63	4.77
	20	110	5.5	10.09
	30	90	3	18.50
	40	70	1.75	31.72
	Number of electrolyte unit (N 1	147	147	0.38
	6	132	22	2.52
	10	120	12	4.62
	18	96	5.33	10.40
	25	75	3	18.50
	32	54	1.68	32.89

This concentration are chosen according to experimental studies where adsorption of electrolyte in zeosils are studying at concentration up to 20 mol/L.

[START_REF] Confalonieri | Structural interpretation of the energetic performances of a pure silica LTA-type zeolite[END_REF][START_REF] Huve | The effect of nanostructures on high pressure intrusion-extrusion of water and electrolyte solutions in hierarchical nanoboxes of silicalite-1[END_REF][START_REF] Ryzhikov | Highpressure intrusion-extrusion of water and electrolyte solutions in pure-silica LTA zeolite[END_REF]

Each system contains 150 particles and the volume of the cubic simulation of side L = 16.646 Å is fixed to match the density of bulk water at atmospheric pressure in absence of ions. The number of electrolyte units N

  TIP4P water molecules (called subsequently water solvation work) while the term (n+m)k B T is the energy cost for removing (n + m) water molecules from vacuum. We previously computed ∆µ CnAm where a single electrolyte unit C n A m is present in the system after water transmutation. However, the solvation free energy of (n+m) water molecules ∆µ (n+m)H 2 O needs to be calculed. Therein, we prepared a box having the same side L as previously, containing 150-(n+m) water molecules and used our modified Widom insertion technique using NCMC moves. We randomly selected (n + m) positions in the simulation box, satisfying the condition that the chosen positions are located at least from a distance of σ Ow /2 from the center of water oxygens in order to avoid initial interpenetration of beads resulting indivergent Lennard-Jones potential and thus blow up of the system during relaxation process by (N, V, E) molecular dynamics. The work required to solvate (n + m) TIP4P water molecules is the excess Widom potential

	defined as:
	∆µ (n+m)H 2 O = -	1 β	ln e -βW (n+m)H 2 O	(150-(n+m))H 2 O

2

O is the work required to solvate (n + m)

Table 2 :

 2 Solvation free energy of electrolytes ∆F CnAm , internal energy ∆U CnAm and entropy ∆S CnAm versus experiments at infinite dilution.

		∆F CnAm (kJ/mol)	∆U CnAm (kJ/mol)	∆S CnAm (J/mol/K)
	C n A m	This work Ref 41	This work Ref 41	This work	Ref 41
	LiCl	-853 ± 3 -849	-886 ± 2 -898	-110 ± 16 -164.1
	NaCl	-738 ± 2 -744	-768 ± 3 -783	-100 ± 14 -132.8
	KCl	-671 ± 2 -671	-695 ± 2 -700	-80 ± 10 -96.7
	RbCl	-658 ± 2 -649	-680 ± 2 -675	-73 ± 10 -85.9
	CsCl	-635 ± 2 -626	-656 ± 2 -650	-70 ± 14 -81.0
	KF	-838 ± 2 -799	-881 ± 2 -847	-143 ± 13 -159.4
	KBr	-658 ± 1 -644	-678 ± 3 -669	-67 ± 12 -80.9
	KI	-603 ± 2 -609	-622 ± 2 -629	-63 ± 12 -62.8
	BaCl 2	-1999 ± 4 -1984	-2061 ± 3 -2065	-206 ± 20 -273.9
	SrCl 2	-2099 ± 3 -2112	-2175 ± 5 -2203	-253 ± 23 -306.3
	CaCl 2	-2217 ± 4 -2239	-2299 ± 3 -2335	-273 ± 23 -324.4
	MgCl 2	-2495 ± 4 -2566	-2582 ± 3 -2682	-290 ± 24 -391.1

  test our implementation of the osmostat in the semi-grand canonical ensemble by imposing a series of chemical potential for CsCl electrolyte ∆µ CsCl to MC simulations starting from bulk water. In SI, supplementary MC simulations were performed by giving extreme values for ∆µ CsCl to probe behaviour of osmostat in asymptotic cases: a highly negative ∆µ CsCl yields a system full of water while setting a highly positive ∆µ CsCl results in a box full of ions. In Figure4, we present the evolution of the concentration of CsCl as a function

	of MC cycles for values of chemical potential of -567.8, -561.2, -551.1 and -525.6 kJ/mol
	(Table S6) that defines a saline reservoir at increasing concentration of respectively 2.41,
	4.77, 10.09 and 18.50 M. Starting from bulk water, the systems in contact with the saline
	reservoir exchange ions and water molecules until equilibrium concentration is reached. A
	higher concentration of the reservoir required more MC steps to drive the system towards
	equilibrium. Only ∼50 MC cycles are required to reach equilibrium in contact with reservoirs
	of 2.41 and 4.77 M while ∼120 and 200 cycles are needed for equilibration with reservoirs of
	10.06 and 18.46 M respectively.

Acknowledgement

We acknowledge funding from Agence Nationale de la Recherche under project MESAMM (ANR-19-CE05-0031).

The choice of P insert and P delete ensures that electrolyte insertion is attempted when no electrolytes remain in the system while electrolyte deletion is attempted when the maximum number of electrolyte units in the system is reached.

• Then, an NCMC move is proposed for alchemically transmute chosen particles during which the alchemical work is recorded in the same fashion as in the section where calculation of chemical potential was performed. Again, the alchemical path involves 200 segments where each segment consists in pertubation through interpolation of nonbonded parameters of the chosen particles, followed by relaxation with 200 NVE MD steps. The alchemical work for transmuting water into an electrolyte unit is the same as defined previously W CnAm while the alchemical work for removing an electrolyte is W -CnAm .

• Let us derive again the probability of accepting and deleting an electrolyte unit C n A m . [START_REF] Ross | Biomolecular simulations under realistic macroscopic salt conditions[END_REF] We recall that the density probability of a system of water and electrolytes in equilibrium with an osmostat at chemical concentration ∆µ CnAm in the semi-grand canonical ensemble writes

The probability of acceptance of a trial NCMC move in which (n+m) water molecules are transmuted into an electrolyte unit, denoted α(N CnAm → N CnAm + 1), is the ratio

Supporting Information Available

Force field parameters, evolution of alchemical work as a function of NCMC moves, error and raw data for alchemical work and chemical potential as well as radial distribution functions, water solvation and tests of the osmostat for extreme cases are provided in supplementary