Maksim Koptelov
email: maksim.koptelov@unicaen.fr

Albrecht Zimmermann
email: albrecht.zimmermann@unicaen.fr

Bruno Crémilleux
email: bruno.cremilleux@unicaen.fr

Link Prediction via Community Detection in Bipartite Multi-Layer Graphs

Keywords: Link prediction, Community detection, Multi-layer graphs

The growing number of multi-relational networks pose new challenges concerning the development of methods for solving classical graph problems in a multi-layer framework, such as link prediction. In this work, we combine an existing bipartite local models method with approaches for link prediction from communities to address the link prediction problem in multi-layer graphs. To this end, we extend existing community detection-based link prediction measures to the bipartite multi-layer network setting. We obtain a new generic framework for link prediction in bipartite multi-layer graphs, which can integrate any community detection approach, is capable of handling an arbitrary number of networks, rather inexpensive (depending on the community detection technique), and able to automatically tune its parameters. We test our framework using two of the most common community detection methods, the Louvain algorithm and spectral partitioning, which can be easily applied to bipartite multi-layer graphs. We evaluate our approach on benchmark data sets for solving a common drug-target interaction prediction task in computational drug design and demonstrate experimentally that our approach is competitive with the state-of-the-art.

Introduction

Many real world applications can be modeled as bipartite graphs, vertices of which are divided into two distinct groups [START_REF] Sun | Neighborhood formation and anomaly detection in bipartite graphs[END_REF]. The problem setting that motivates our work is the prediction of links between drug candidates and biological targets, an essential step of computational drug design. But there are other link prediction settings that fall into the same category, for instance user-product recommendation.

The available data on drug-target interaction prediction are of heterogeneous structure, i.e. represented by networks the edges of which have different origins, thus making the use of most existing link prediction methods in a straightforward way impossible. Current solutions are limited by the number or type of networks, often referred to as layers, e.g. three layers with two assumed to be similarity networks [START_REF] Chen | Drug-target interaction prediction by random walk on the heterogeneous network[END_REF][START_REF] Lim | Improved genome-scale multi-target virtual screening via a novel collaborative filtering approach to cold-start problem[END_REF][START_REF] Buza | Aladin: A new approach for drug-target interaction prediction[END_REF].

To address these restrictions, we take inspiration from existing methods that use community detection to perform link prediction [START_REF] Guimerà | Missing and spurious interactions and the reconstruction of complex networks[END_REF][START_REF] Yan | Finding missing edges in networks based on their community structure[END_REF][START_REF] Soundarajan | Using community information to improve the precision of link prediction methods[END_REF][START_REF] Valverde-Rebaza | Link prediction in online social networks using group information[END_REF]. While this decouples the problem into how to find communities in multilayer graphs, and how to exploit them for prediction, existing link prediction measures [START_REF] Cannistraci | From link-prediction in brain connectomes and protein interactomes to the local-community-paradigm in complex networks[END_REF][START_REF] Xie | Potential links by neighbor communities[END_REF][START_REF] Ding | Prediction of missing links based on community relevance and ruler inference[END_REF] are not directly applicable to the bipartite setting. To address this restriction, we extend several of those measures to our problem setting. In addition, we go a step further by proposing alternatives to those measures based on an existing bipartite local model. While we evaluate two concrete approaches for community detection, spectral partitioning and the Louvain algorithm, both of which can be easily applied to multi-layer graphs, we do not require any particular community detection approach. Our long-term contribution is the adaptation of existing link-prediction-by-community-detection measures to the bipartite multilayer setting, which we evaluate experimentally, and the selection of best measure and community detection approach combination. In addition, we demonstrate that the parameter settings of the spectral partitioning method can be effectively set via internal cross-validation.

The rest of the paper is organized as follows. Section 2 provides basic notations and definitions. Section 3 discusses related work on link prediction and community detection in multi-layer networks. Section 4 explains how we adapt existing measures for our framework. Section 5 describes the data used for evaluation, the experimental setup and presents the results. Finally, Section 6 concludes and outlines the future work.

Definitions

A graph is a tuple G =< V, E >, where V = {v 1 , v 2 , ..., v n } denotes a set of vertices or nodes, and E ⊆ V × V a set of edges defined by distinct vertex pairs (u, v) ∈ V × V with u = v (without self-loops). We also use the notion of a bipartite graph, which we define as a graph the vertices of which can be divided into two classes V 1 and V 2 such that there is no edge between vertexes of the same class:

G =< V 1 ∪ V 2 , E >, E ⊂ V 1 × V 2 .
We address weighted and unweighted graphs in the same manner. We define a weighted graph as one with a labeling function for edges E → A e with A e ∈ [0, 1], where 0 means no interaction between vertices, 1 confirmed interaction, and an intermediate value represents interaction probability. An unweighted graph is one where every edge is labeled by 1.

To exploit different sources of information in one single structure, we employ multi-layer networks. We define a multi-layer network as a weighted graph where more than one edge (u, v) can exist for a pair of vertices u, v. Multi-layer networks can be decomposed into disjunct set of graphs G l that contain at most a single edge for each pair of vertices, called layers or just networks. As we wrote above, our original setting is a bipartite one. To combine it with the multilayer framework, we define a bipartite multi-layer graph as a multi-layer network the vertices of which can be divided into two classes, and where exactly one of the layers is a bipartite graph. Note, according to the classification of Kivel et al. [START_REF] Kivelä | Multilayer networks[END_REF], the multi-layer networks used in this work are not node-aligned 1 , not layer-disjoint2 , have diagonal couplings3 which are categorical4 , and the number of layers can be any.

We represent graphs as matrices. The adjacency matrix A has size n × n, n = |V |, and A ij represents the weight of the edge (v i , v j ). In the case of multilayer networks, we aggregate the weights of multiple edges between v i and v j by summing them up. Note, A has zeros on the main diagonal, because graphs as used in this work have no self-loops. The degree matrix D is the diagonal matrix

D =    deg(v 1 ) 0 . . . 0 . . . 0 0 deg(v n )    , deg(v i ) = n j=1
w ij of same size as A, where deg(v i ) represents the degree of vertex v i . The degree of a vertex is the sum of the weights of the edges adjacent to v i [START_REF] Gallier | Notes on elementary spectral graph theory. applications to graph clustering using normalized cuts[END_REF].

The last, and arguably most important, matrix used in this paper is the Laplacian matrix. The Laplacian matrix, denoted by L, is a matrix of the same dimensionality as A and D, defined as the difference between the degree matrix and the adjacency matrix: L = D -A. L has the same values as D on the diagonal, and off the diagonal L ij is equal to -A ij .

Related work

Existing methods for link prediction in bipartite multi-layer networks for addressing the drug-target interaction problem can be grouped into three classes: similarity based, random-walker based, and latent models based. The first group assumes 2 out of 3 possible layers to be similarity networks for drugs and targets respectively, and exploits similarity information to perform link prediction on the third bipartite layer [START_REF] Ding | Similarity-based machine learning methods for predicting drug-target interactions: a brief review[END_REF][START_REF] Buza | Aladin: A new approach for drug-target interaction prediction[END_REF]. The second models the behavior of a random-walker to perform link prediction in multi-layer graph using PageRank adaptations [START_REF] Chen | Drug-target interaction prediction by random walk on the heterogeneous network[END_REF][START_REF] Cheng | Prediction of chemical-protein interactions network with weighted network-based inference method[END_REF][START_REF] Koptelov | Link prediction in multi-layer networks and its application to drug design[END_REF]. Such methods are dependent on fixing the similarity networks and while we extended the approach to any number of networks in [START_REF] Koptelov | Link prediction in multi-layer networks and its application to drug design[END_REF], it pays for this flexibility with high computational cost. The last group of methods maps drugs, targets and their interactions into a combined feature space, and performs drugtarget interaction prediction using distance functions or regression analysis [START_REF] Yamanishi | Prediction of drug-target interaction networks from the integration of chemical and genomic spaces[END_REF][START_REF] Zheng | Collaborative matrix factorization with multiple similarities for predicting drug-target interactions[END_REF]. The most recent family of methods in this mold is often referred to as graph embeddings [START_REF] Goyal | Graph embedding techniques, applications, and performance: A survey[END_REF]. The main disadvantage of this group of methods is a certain lack of interpretability.

The idea to use community information to predict links in graphs is not novel. Clauset et al. [START_REF] Clauset | Hierarchical structure and the prediction of missing links in networks[END_REF] proposed to exploit a learned hierarchical generative community model to estimate the probabilities of missing links in partially known networks.

The authors of [START_REF] Guimerà | Missing and spurious interactions and the reconstruction of complex networks[END_REF][START_REF] Yan | Finding missing edges in networks based on their community structure[END_REF][START_REF] Valverde-Rebaza | Link prediction in online social networks using group information[END_REF] combine community detection with existing edge prediction methods to improve the prediction accuracy. These methods are based on the hypothesis that vertices in the same community have similar properties, and missing edges are more likely to be found within communities than elsewhere. Missing edges are predicted by node similarity using nearest-neighbor measures [START_REF] Yan | Finding missing edges in networks based on their community structure[END_REF], Stochastic Block Models [START_REF] Guimerà | Missing and spurious interactions and the reconstruction of complex networks[END_REF], or in-group and out-group neighbor similarity measures [START_REF] Valverde-Rebaza | Link prediction in online social networks using group information[END_REF]. Edges can be predicted for vertices belonging to the same community even if there is no path between them within the community [START_REF] Jalili | Link prediction in multiplex online social networks[END_REF]. The density of links in a particular community or between two communities can be exploited in a naïve Bayes model to predict links [START_REF] Liu | Correlations between community structure and link formation in complex networks[END_REF]. The authors of [START_REF] Ahn | Link communities reveal multiscale complexity in networks[END_REF] reimagine communities as groups of edges rather than vertices, and [START_REF] Soundarajan | Using community information to improve the precision of link prediction methods[END_REF] use community detection to modify similarity measures. Finally, there is a set of methods which extend the concept of shared neighborhoods [START_REF] Liben-Nowell | The link-prediction problem for social networks[END_REF] to community neighborhoods [START_REF] Cannistraci | From link-prediction in brain connectomes and protein interactomes to the local-community-paradigm in complex networks[END_REF][START_REF] Xie | Potential links by neighbor communities[END_REF][START_REF] Ding | Prediction of missing links based on community relevance and ruler inference[END_REF]. In addition, in [START_REF] Hristova | A multilayer approach to multiplexity and link prediction in online geo-social networks[END_REF] neighborhood measures have been extended to multiple layers.

Community detection in multi-graphs can also be performed in different ways: directly, by ensemble-based methods, or by graph flattening. The direct methods perform discovery of communities on the multi-layer network directly, e.g. by adapting objective functions for community detection to the multi-layer setting [START_REF] Kuncheva | Community detection in multiplex networks using locally adaptive random walks[END_REF][START_REF] Tagarelli | Ensemble-based community detection in multilayer networks[END_REF][START_REF] De Bacco | Community detection, link prediction, and layer interdependence in multilayer networks[END_REF]. Ensemble-based methods perform community detection on each layer separately, and aggregate discovered communities afterwards [START_REF] Tagarelli | Ensemble-based community detection in multilayer networks[END_REF]. Flattening approaches, finally, summarize multiple edges into single ones and use the resulting single-layer network to discover communities by using one of the common community detection approaches such as spectral partitioning [START_REF] Leskovec | Mining of massive datasets[END_REF] or Louvain algorithm [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF]. In this work we use the last type of approaches due to their ease of use and potentially low computational complexity.

Our approach

In our problem setting, we want to predict links between two distinct types of nodes, e.g. drugs and targets in our experiments. To achieve this, we perform communities discovery using en exsting community detection approach, then exploit the discovered communities to solve the link prediction task.

Link prediction by community detection in bipartite setting

To be able to use existing link-prediction-by-community-detection measures we have to adapt them to the bipartite setting. Due to the construction of the networks we use and the community detection methods we evaluate, resulting communities can be mixed, i.e. containing both types of nodes, drugs and targets, as well as pure, of either type, drugs or targets only. Also, the community detection methods we use produce non-overlapping communities only. Mixed communities can be exploited directly with existing measures for link prediction via community detection (see Section 3), but pure communities cannot, ignoring a large number of drug-target pairs. To overcome this, we treat all communities as non-mixed and split mixed communities into pure ones. Notably, this split does not have to be done explicitly, but a mixed community can be treated as two pure ones with links between them. We exploit discovered communities in one of the two proposed ways: by matching "community to community" or "node to community".

Community to community In this case, each drug community is paired with each target community, then an adapted measure used to perform link prediction between paired communities. Each non-interacting drug-target pair between paired communities is assigned the same link probability score. At the end of the matching, each non-interacting drug-target pair from the network will have been assigned a single score, which can be used to rank predictions. We refer to this approach as community to community (or CC ) formulation.

Node to community Another way of exploiting communities is to pair each node of one type with communities of the other type. The advantage of that method is that for a selected drug d i and target t j , the prediction can be made twice: once analyzing connections of a drug with target communities and second analyzing target connections with drug communities, providing a more reliable estimate. The approach based on this idea is called Bipartite Local Models [START_REF] Bleakley | Supervised prediction of drug-target interactions using bipartite local models[END_REF]. The link probability score between d i and t j is computed by aggregating the two results [START_REF] Buza | Aladin: A new approach for drug-target interaction prediction[END_REF]. We report results using mean as an aggregation function. Our experiments showed that the difference between max and mean is negligible, and we use mean to get a more reliable result. We do not consider min as aggregator, because in case of no evidence for existence of the link in one of the independent predictions the combined probability is also 0. We refer to this approach as node to community (or NC ) formulation.

Existing link prediction measures adaptation

We divide all existing link prediction measures into two categories: neighborhood measures and others, which we refer as community-based. The first group of measures are based on the notion of neighborhood, i.e. the set of vertices directly connected to the examined vertices. The semantic similarity between neighborhood and community, i.e. sets of vertices in both cases, allows us to use neighborhood measures in our setting. The other measures are not based on a notion of neighborhood, but on other metrics, and thus are grouped into separate group in our work.

Neighborhood measures Many existing link prediction measures exploit the neighborhoods of vertices e.g. in the form of common neighbors (CN), the Jaccard coefficient (JC), a preferential attachment measure (PA), or SimRank (SR) [START_REF] Liben-Nowell | The link-prediction problem for social networks[END_REF]:

CN (d i , t j ) = |{v | (d i , v) ∈ E} ∩ {u | (t j , u) ∈ E}|, (1)
JC(d i , t j ) = CN (d i , t j ) |{v | (d i , v) ∈ E} ∪ {u | (t j , u) ∈ E}| , (2)
P A(d i , t j ) = |Γ (d i )| • |Γ (t j )|, |Γ (v)| = deg(v), (3)
SR(d i , t j ) = CN (d i , t j ) P A(d i , t j ) . (4)
Due to the nature of communities we obtain, there is little overlap between vertices' neighborhoods, preventing the direct use of neighborhood-based measures.

To overcome this, we adapt neighborhood measures for use with our communities, treating them like neighborhoods: the CN measure turns the number of common neighbors of communities d i and t j into the number of connections, JC represents the fraction of all possible connections of d i and t j that are connected to both, PA is defined by a product of degrees of communities d i and t j , finally SR is equal to the number of connections of communities d i and t j normalized by the product of their degrees. Using the CC and NC formulations our bipartite adaptations take the form:

1. Instead of the measures from Eq. 1-4 we define CN CC , JC CC , PA CC and SR CC versions corresponding to CC matching:

CN CC (d i , t j ) = |{(d, t) ∈ E | d ∈ C(d i ), t ∈ C(t j )}|, (5)
JC CC (d i , t j ) = CN CC (d i , t j ) |C(d i )| • |C(t j )| . (6)
P A CC (d i , t j ) = |Γ C (d i )| • |Γ C (t j )|, (7)
SR CC (d i , t j ) = CN CC (d i , t j ) P A CC (d i , t j ) , (8)
where Γ C represents the neighborhood of a community, the degree of which we defined in accordance with [START_REF] Clauset | Hierarchical structure and the prediction of missing links in networks[END_REF] as

|Γ C (v)| = d l ∈C(v),d k / ∈C(v) w(d l , d k ).
2. The NC version of CN is defined as the average of the two independent predictions,

CN N C (d i ) = |{t | (d i , t ∈ E, t ∈ C(t j )}|, CN N C (t j ) = |{d | (t j , d) ∈ E, d ∈ C(d i )
}|, for d i , and t j respectively:

CN N C (d i , t j ) = 1 2 (CN N C (d i ) + CN N C (t j )) , (9)
In the same manner, we can define NC versions for other measures from [START_REF] Liben-Nowell | The link-prediction problem for social networks[END_REF], the Jaccard coefficient (Eq. 10), preferential attachment (Eq. 11) and SimRank (Eq. 12), taking into account that

P A N C (d i ) = d k A(d i , d k ) • |Γ C (t j )| and P A N C (t j ) = t k A(t j , t k ) • |Γ C (d i )|: JC N C (d i , t j ) = 1 2 CN N C (d i ) |C(t j )| + CN N C (t j ) |C(d i )| . ( 10
) P A N C (d i , t j ) = 1 2 (P A N C (d i ) + P A N C (t j )), (11)
SR N C (d i , t j ) = 1 2 CN N C (d i ) P A N C (d i ) + CN N C (t j ) P A N C (t j ) , (12)
Community-based measures Other measures proposed in the literature are based on one or several of the following assumptions: all vertices have the same semantic, all edges have the same semantic, edges are unweighted, or vertices whose link is to be predicted find themselves in the same community. We therefore cannot use most of the measures proposed in the literature but we can adapt some to our bipartite setting.

1. Cannistraci et al. [START_REF] Cannistraci | From link-prediction in brain connectomes and protein interactomes to the local-community-paradigm in complex networks[END_REF] in their CAR-based measures propose to exploit the density of communities to reward (or penalize) densely (sparsely) connected neighbors of the vertices whose link is to be predicted. Our adapted CARbased common neighbors (CCN) will be defined as CN regularized by community local degree, which is in turn defined as the sum of weights of all edges inside community. The NC formulation of CCN takes a form:

CCN N C (d i , t j ) = 1 2 (CCN N C (d i ) + CCN N C (t j )) , (13)
with CCN N C (d i ) and CCN N C (t j ) in turn defined as:

CCN N C (d i ) = |{t | (d i , t) ∈ E, t ∈ C(t j )}| • t l ,t k ∈C(tj ) A(t l , t k ), CCN N C (t j ) = |{t | (t j , d) ∈ E, d ∈ C(d i )}| • d l ,d k ∈C(di) A(d l , d k ).
In the same manner, the CAR-based Jaccard coefficient (CJC) is redefined as CCN normalized by the size of the community:

CJC N C (d i , t j ) = 1 2 (CJC N C (d i ) + CJC N C (t j )) , (14)
CJC N C (d i ) = CCN N C (d i ) |C(t j )| , CJC N C (t j ) = CCN N C (t j ) |C(d i )| .
2. Xie et al. [START_REF] Xie | Potential links by neighbor communities[END_REF] propose to exploit the connection of vertices to communities, summing over all communities. Our adaptation of their measure, which we refer to as Neighboring community-based (NCB) is defined as the normalized sum of all CN regularized by the size of the respective community. The NC formulation of this measure will take the form:

N CB N C (d i , t j ) = 1 2 (N CB N C (d i ) + N CB N C (t j )) , (15)
with N CB N C (d i ) and N CB N C (t j ) in turn:

N CB N C (d i ) = ct k=1 |{t | (d i , t) ∈ E, t ∈ C k }| |C k | • |{t | t ∈ C k }| |T | , N CB N C (t j ) = c d k=1 |{d | (t j , d) ∈ E, d ∈ C k }| |C k | • |{d | d ∈ C k }| |D| .
Moreover, assuming communities are pure, i.e. consisting only of either drugs or targets, these equations can be simplified to the sum of all CN normalized by the number of vertices of one type:

N CB N C (d i ) = 1 |T | ct k=1 |{t | (d i , t) ∈ E, t ∈ C k }|, N CB N C (t j ) = 1 |D| c d k=1 |{d | (t j , d) ∈ E, d ∈ C k }|.
3. Ding et al. [START_REF] Ding | Prediction of missing links based on community relevance and ruler inference[END_REF], finally, propose to exploit the neighborhoods of communities.

Our adaptation of their measure, which the authors refer to as Community relevance Jaccard coefficient (CRJC), is defined as a number of common nodes of examined communities and nodes of the opposite type connected to those communities normalized by the total number of nodes in this selection. The adapted measure is better suited to a CC formulation:

CRJC CC (d i , t j ) = |CRJC CC (d i ) ∩ CRJC CC (t j )| |CRJC CC (d i ) ∪ CRJC CC (t j )| , (16)
with CRJC CC (d i ) and CRJC CC (t j ):

CRJC CC (d i ) = {t | (d, t) ∈ E, d ∈ C(d i )} ∪ {d | d ∈ C(d i )}, CRJC CC (t j ) = {d | (t, d) ∈ E, t ∈ C(t j )} ∪ {t | t ∈ C(t j )}.

Experimental evaluation

To evaluate the two community detection algorithms, the effect of their parameter settings, and the prediction measures defined in the preceding section, we performed experiments on a number of benchmark data sets for ligand-target activity prediction.

Experimental protocol

We begin by evaluating the different measures described in Section 4 with two common community detection approaches, keeping most of the parameters fixed, and select the best performing measure. Following this, we show how an internal cross-validation can be used to fixed the methods' parameters, and report on their results. Finally, we show what happens if we optimized those parameters on the test data, giving us the arguably best results, and show that the internal cross-validation creates models of similar quality, i.e. using cross-validation is an appropriate method to fix parameters.

Community detection methods

We test our approach with spectral partitioning [START_REF] Leskovec | Mining of massive datasets[END_REF] and the Louvain algorithm [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF] as community detection approaches. The first finds the best cut to partition nodes based on eigenvalues of the Laplacian matrix and a threshold method [START_REF] Spielman | Spectral partitioning works: Planar graphs and finite element meshes[END_REF][START_REF] Fortunato | Community detection in graphs[END_REF], while the second greedily optimizes modularity, generic measure to determine the quality of each partition produced by a community detection approach [START_REF] Newman | Finding and evaluating community structure in networks[END_REF][START_REF] Fortunato | Community detection in graphs[END_REF]. We apply spectral partitioning to multi-layer graphs by flattening the graph, i.e. summing edge weights to derive the adjacency and degree matrices before performing partitioning. We apply the Louvain algorithm to multi-layer graphs by flattening the graph as well. Since the algorithm is not limited to adjacency and degree matrices, we create an instance of a single graph, multiple edges of which are aggregated with the sum function.

Parameters to optimize Spectral partitioning has two parameters: m value and the thresholding method. The m parameter represents the number of eigenvectors corresponding to the m smallest non-zero eigenvalues used to partition the graph into at most 2 m groups. As thresholds methods we can use sign cut or default, which partitions entries based on whether they are greater or less than zero, bisection cut or median, using the median value of entries in an eigenvector as a threshold, producing two components of approximately equal size [START_REF] Guattery | On the performance of spectral graph partitioning methods[END_REF]. We also evaluate mean, which uses the average, and sum that exploits the fact that there are approximately equally as many positive and negative values in eigenvectors of Laplacian matrix, such that in practice they sum to a value close to zero. Moreover, the same threshold can be applied to all eigenvectors, or each individual eigenvector can have its own threshold. We call the former approach global, and the latter individualized. Additionally, the global threshold can be computed by applying the aggregating function (mean, median or sum) to all eigenvectors or only to the m actually used. We refer to this latter type as localized. To sum up, we evaluate 9 different thresholding methods: global, localized, individualized and their combinations with mean, median and sum. The combination global sum is a special case since taking the first eigenvector, whose entries all have the same, positive, value, into account violates the close to zero property sum thresholding exploits. We therefore do not evaluate that thresholding method, but add default thresholding to the mix for the experimental evaluation.

The Louvain algorithm has only one parameter -resolution limit -which defines a modularity scale. Practically speaking, at different moments of time t, the difference between optimal partitioning and partitioning produced by the Louvain is different and the resolution parameter represents this change in time [START_REF] Lambiotte | Laplacian dynamics and multiscale modular structure in networks[END_REF].

Data sets We perform our experiments on the data sets introduced in [START_REF] Yamanishi | Prediction of drug-target interaction networks from the integration of chemical and genomic spaces[END_REF]: Enzyme, G-protein coupled receptors (GPCR), Ion Channels (IC) and Nuclear Receptors (NR). In addition, we use the Kinase set [START_REF] Davis | Comprehensive analysis of kinase inhibitor selectivity[END_REF]. These data sets have been used in prior work on drug-target interaction prediction [START_REF] Chen | Drug-target interaction prediction by random walk on the heterogeneous network[END_REF][START_REF] Zheng | Collaborative matrix factorization with multiple similarities for predicting drug-target interactions[END_REF][START_REF] Lim | Improved genome-scale multi-target virtual screening via a novel collaborative filtering approach to cold-start problem[END_REF][START_REF] Buza | Aladin: A new approach for drug-target interaction prediction[END_REF], and can be considered benchmarks. The data consist of 3 networks: drug similarities, target similarities and drug-target interaction (the bipartite graph). The interaction network is presented as binary relation lists, while the similarity data are presented by matrices. The data were transformed to adjacency lists, and for the sake of (computational) convenience we mapped ligand and drug names to consecutive integers, i.e. l0, l1, l2, . . ., t0, t1, t2, . . ., the required input for our implementation. The data sets' basic properties are presented in Table 1. Evaluation protocol To perform evaluation of our experiments we performed a 5×5-fold cross-validation, with each fold containing 20% of all drug-target interactions, acting as test set for link prediction once, while community detection is performed on the other 80%. The process is repeated 5 times, the results are averaged among all runs.

Quality measures

We evaluate all the predictions by AUC (Area Under ROC Curve) and AUPR (Area Under Precision-Recall Curve), averaging the results. We also report standard deviation values when it is applicable.

Implementation We implemented spectral partitioning and all link prediction measures in Python 5 , using the networkx library to model the multi-layer network, the python-louvain package as Louvain algorithm implementation, NumPy for all matrix computations and sklearn to calculate the curves.

Experimental results

Link prediction measures evaluation We first test the different link prediction measures from Section 4 on communities produced by either spectral partitioning or the Louvain algorithm. To reduce computational complexity, we use default parameters of the community detection approaches: default as threshold for spectral partitioning and 1.0 as resolution limit for the Louvain algorithm. Note, that we optimize m for spectral partitioning on the test data in this experiment, because there is no a priori number of eigenvalues that will fit all data. The results are presented in Fig. 1, with community to community formulations on the left, node to community ones on the right of each plot. The best-performing measure for each group is indicated by a + sign over the corresponding bar. We also report on the figure the optimal m value for every measure in spectral partitioning.

The results show that a number of measures, e.g. SR CC , CN N C , SR N C , CCN N C , CJC N C , N CB N C , have acceptable performance in terms of AUC on most of data sets while the Jaccard coefficient usually performs best for both the CC and NC versions. These two, JC CC and JC N C , are therefore the measures we will use going forward. Another result is that for these parameter settings spectral partitioning and the Louvain algorithm give approximately the same AUC, but the latter improves on AUPR. Finally, using node to community predictions requires a lower m, i.e. less fine-grained partitions, for spectral partitioning.

Parameter selection for spectral partitioning via internal cross-validation In link prediction, as in any other predictive task, the main issue is choosing parameter values, in the case of spectral partition m and the thresholding method. In the absence of other knowledge, one would use cross-validation as a systematic method to optimize parameters. Several-fold cross-validation is also the method of choice to evaluate the performance of a classifier, however, so that we use a double cross-validation in this section: splitting off an external test fold (containing 20% of present edges) to evaluate the model, and using an internal five-fold cross-validation to fix the model's parameters.

Table 2 presents the results. It shows both the results of internal evaluation, i.e. on the validation set used to fix parameter values, and of the external evaluation, i.e. on the unseen training data. The main conclusion to draw is that those values align very closely, i.e. that there is not risk of overfitting when building the model.

Baseline comparison

We have addressed this problem setting in prior work, using a random walk approach [START_REF] Koptelov | Link prediction in multi-layer networks and its application to drug design[END_REF] which is basically an extension of PageRank for any number of layers, and thus can be considered as a baseline in this work. JC CC and JC N C in combination with both spectral partitioning and the Louvain clearly outperform that baseline in terms of AUC (0.84, 0.8, 0.76, 0.63 and 0.61 for the Enzyme, GPCR, IC, NR, and Kinase data sets respectively) and AUPR (0.15, 0.22, 0.27, 0.26 and 0.13 for the Enzyme, GPCR, IC, NR, and Kinase data sets respectively).

Performance ceiling and state of the art comparison While we have shown above that using an internal cross-validation can be used to use our method's parameters, one could ask the question whether those results approach the best possible results on the data, and how they compare to the state of the art. To explore this issue, we performed parameter optimization on the test data, an unrealistic setting since we take the label information of unseen data into account. The performance in this setting can be interpreted as the ceiling of what is achievable. Tables 3 and4 list optimal parameter choices for the two community detection methods, as well as predictive performance. The results of Table 3 are very close to those of Table 2, showing that the internal cross validation gets close to the achievable ceiling. 6 We also report standard deviations, which show the model's performance is very stable. As one can see, the NC setting continues to outperform CC for all data sets. In addition, the Louvain algorithm provides not only better AUC than spectral partitioning, but significantly better AUPR, different from the results for the default parameter setting (Fig. 1).

Concerning comparison with the state-of-the-art, the Louvain using JC N C comes close to the performance of the state-of-the-art methods reported in [START_REF] Buza | Aladin: A new approach for drug-target interaction prediction[END_REF] in terms of AUC: 0.97, 0.95, 0.98, 0.88 and 0.9 for the Enzyme, GPCR, IC, NR, and Kinase data sets respectively. Interpretability Our approach offers a straightforward option for interpretation/explanainability of a link prediction: for each of the two vertices, we can show the communities they belong to, the weights of intra-community edges, the number and layout of inter-community edges, and their numerical translation by the measure and matching technique. This is a possibility that is not available for recent, well-performing techniques based on graph embeddings.

Scalability In order to verify scalability of our approach we tested it on a bigger data set, IUPHAR, having 6 layers, and described in detail in [START_REF] Koptelov | Link prediction in multi-layer networks and its application to drug design[END_REF]. The data set basic properties are presented in Table 1. We used spectral partitioning with default threshold and an optimal value of m = 400 optimized on the test data. We used JC, the best performing measure and CC as a matching technique for performance reasons. As for the smaller networks, the prediction quality (AUC=0.7365 ± 0.02 std, AUPR 0.01 ± 0.0 std) is much better than in our previous work [START_REF] Koptelov | Link prediction in multi-layer networks and its application to drug design[END_REF] (AUC 0.5735) (which was derived from leave-out-out crossvalidation). Running times for a single test fold, the majority of which is taken up by community matching, are shown in Table 1. We can see that the quotient of the number of IUPHAR edges to the number in smaller networks is lower than that of running times, indicating that the method scales when using spectral partitioning.

Conclusion and perspectives

We have presented a framework for link prediction in bipartite multi-layer graphs using graph community structure and link prediction measures adapted from those proposed in the literature. We have found empirically that combining the well-known and relatively straightforward Jaccard coefficient, particularly in a BLM formulation, with the Louvain algorithm for community detection allows us to achieve results that are competitive with the state-of-the-art.

We have limited ourselves to two easy-to-use community detection methods in this work, and will evaluate the use of other methods in the future. We also intend to perform experiments on larger data sets that have shown themselves to be too computationally expensive for methods such as the one proposed in [START_REF] Buza | Aladin: A new approach for drug-target interaction prediction[END_REF]. Finally, we intend to add layers derived from other information sources to the networks and use our approach to identify possible redundancies among them.

Fig. 1 .

 1 Fig. 1. Link prediction measures evaluation on the benchmark data sets. The symbol + denotes the best performing measure for each group of formulations in each data set.

Table 1 .

 1 Basic properties of Benchmark and IUPHAR data sets and running times (* -for 1 fold in average with Spectral partitioning and JCCC measure)

	Data set Drugs Targets Interactions Layers |V |	|E|	Sparsity CC Running
										time*, s
	Enzyme 445	664	2926	3	1109 321832	0.524	1	58.67
	GPCR 223	95	635	3	318	29853	0.592	1	3.99
	IC	210	204	1476	3	414	44127	0.516	1	7.06
	NR	54	26	90	3	80	1846	0.584	1	0.15
	Kinase	68	442	1527	3	510 101266	0.780	1	9.94
	IUPHAR 8137 2502	12456	6	10639 26706838 0.472	1 3477.8

Table 2 .

 2 Spectral partitioning parameters optimization via Internal cross-validation

	Data set Measure	Internal validation	External validation
			µ AUC σ µ AUPR σ µ AUC σ µ AUPR σ
	Enzyme JCCC	0.85 0.00 0.14 0.04 0.85 0.01 0.19 0.07
		JCNC	0.91 0.00 0.19 0.02 0.92 0.01 0.26 0.06
	GPCR JCCC	0.79 0.01 0.11 0.02 0.80 0.02 0.15 0.02
		JCNC	0.84 0.01 0.19 0.01 0.85 0.01 0.25 0.02
	IC	JCCC	0.82 0.01 0.31 0.03 0.83 0.02 0.40 0.04
		JCNC	0.87 0.00 0.32 0.06 0.88 0.01 0.41 0.05
	NR	JCCC	0.73 0.02 0.21 0.03 0.72 0.06 0.24 0.08
		JCNC	0.75 0.03 0.19 0.02 0.77 0.06 0.23 0.13
	Kinase JCCC	0.76 0.01 0.17 0.01 0.77 0.03 0.23 0.02
		JCNC	0.85 0.00 0.26 0.02 0.86 0.01 0.35 0.02

Table 3 .

 3 Spectral partitioning threshold and m value optimization

	Data set Measure	Optimal parameters	Performance
			threshold	m value AUC σ AUPR σ
	Enzyme JCCC localized median	11	0.86 0.02 0.15 0.04
		JCNC individual mean	5	0.92 0.01 0.29 0.05
	GPCR JCCC localized median	17	0.80 0.03 0.17 0.03
		JCNC	global mean	3	0.85 0.02 0.27 0.04
	IC	JCCC individual median	13	0.84 0.03 0.47 0.09
		JCNC localized median	7	0.89 0.02 0.43 0.06
	NR	JCCC	global mean	4	0.77 0.08 0.29 0.11
		JCNC	individual sum	2	0.78 0.05 0.25 0.09
	Kinase JCCC localized median	11	0.78 0.02 0.24 0.03
		JCNC	localized mean	6	0.86 0.01 0.36 0.03

Table 4 .

 4 Louvain algorithm resolution optimization

	Data set Measure Optimal	Performance
			resolution AUC σ AUPR σ
	Enzyme JCCC	0.4	0.95 0.01 0.68 0.03
		JCNC	0.6	0.96 0.01 0.71 0.02
	GPCR JCCC	0.7	0.83 0.03 0.35 0.05
		JCNC	0.8	0.89 0.02 0.47 0.08
	IC	JCCC	0.4	0.95 0.01 0.68 0.03
		JCNC	0.8	0.97 0	0.79 0.03
	NR	JCCC	0.8	0.83 0.06 0.31 0.09
		JCNC	0.9	0.84 0.06 0.37 0.09
	Kinase JCCC	0.7	0.73 0.01 0.19 0.01
		JCNC	0.9	0.91 0.01 0.5 0.03

All nodes are shared between all layers

Each node is present only in a single layer

Inter-layer edges, that cross layers, are only between nodes and their counterparts

Diagonal couplings for which all possible inter-layer edges are present

We had to rerun internal parameter optimization for Louvain and the results were not yet available at the time of writing.


A J L Garcia

J.-F Gonzalez

Evolution of porous dust grains in protoplanetary discs -I. Growing grains

Keywords: Protoplanetary discs -Hydrodynamics -Planets and satellites, formation -Methods, numerical

One of the main problems in planet formation, hampering the growth of small dust to planetesimals, is the so-called radial-drift barrier. Pebbles of cm to dm sizes are thought to drift radially across protoplanetary discs faster than they can grow to larger sizes, and thus to be lost to the star. To overcome this barrier, drift has to be slowed down or stopped, or growth needs to be sped up. In this paper, we investigate the role of porosity on both drift and growth. We have developed a model for porosity evolution during grain growth and applied it to numerical simulations of protoplanetary discs. We find that growth is faster for porous grains, enabling them to transition to the Stokes drag regime, decouple from the gas, and survive the radial-drift barrier. Direct formation of small planetesimals from porous dust is possible over large areas of the disc.

INTRODUCTION

In protoplanetary discs, planetesimals are thought to be the building blocks of planets [START_REF] Safronov | Evolution of the protoplanetary cloud and formation of the earth and planets[END_REF]. However, our understanding of the collisional growth of sub-micron-sized monomers up to kilometre-sized planetesimals is hindered by problems called 'barriers' in the planet formation theory. Because of the gas pressure gradient, a differential velocity exists between gas and dust, and grains experience an aerodynamic drag force. Grains thus lose angular momentum, making them settle down to the mid-plane and drift inwards. The influence of the drag force can be measured by the Stokes number St, i.e. the ratio between the grain stopping time (the time for a grain to reach the gas velocity) to the Keplerian orbital time. Small grains, with St 1, and large grains, with St 1, are respectively very coupled with the gas and weakly affected by the gas drag and so, drift slowly. Grains with intermediate sizes, i.e. St∼1, are subject to the fastest radial drift. Thus, initially small grains growing to this intermediate regime can be accreted onto the star faster than they can grow to reach planetesimal sizes (Adachi et al. 1976;Weidenschilling 1977). To overcome this barrier, called the "radial-drift barrier", grain growth needs to be accelerated or drift slowed down. So-called dust traps have been shown to be a good way to slow down drift by accumulating grains in a gas pressure maximum (Haghighipour 2005) such as, e.g., in vortices (Barge & Sommeria 1995;Meheut et al. 2012;Zhu et al. 2014), snow lines (Kretke & Lin 2007;Brauer et al. 2008; Dra ¸żkowska E-mail: anthony.garcia@ens-lyon.fr † E-mail: jean-francois.gonzalez@ens-lyon.fr & Alibert 2017) or planet gaps (Paardekooper & Mellema 2004;Fouchet et al. 2007Fouchet et al. , 2010;;Ayliffe et al. 2012;Pinilla et al. 2012;Zhu et al. 2012Zhu et al. , 2014)). More recently, Gonzalez et al. (2017a,b) showed that self-induced dust traps can form without an initially present pressure maximum in the disc, where growing and fragmenting grains can accumulate thanks to the combined effect of the back-reaction of dust on gas and large-scale gradients.

In this work, we choose to investigate on the growth of porous grains as a solution. Indeed, at a given mass, porous grains have a larger collisional cross-section than compact dust, which can accelerate growth. Until now, grains have been considered as compact in most studies of grain growth (e.g. Brauer et al. 2008;Laibe et al. 2008;Dr ążkowska et al. 2014;Gonzalez et al. 2015a), for the sake of simplicity and because of a lack of knowledge on how porosity evolves with collisions. However, some objects in the Solar system such as comets have appeared to be porous (Blum et al. 2006;A'Hearn 2011;Sierks et al. 2015). Experiments have also shown that low-velocity collisions lead to the formation of porous aggregates (Blum 2004). Furthermore, experimental and numerical works have studied faster collisions and found that collisional energy is dissipated by internal restructuring and that compression can occur (Dominik & Tielens 1997;Blum & Wurm 2000;Wada et al. 2007Wada et al. , 2008Wada et al. , 2009;;Suyama et al. 2008Suyama et al. , 2012;;Seizinger et al. 2012). The importance of porosity on the collisional evolution of dust grains, and in particular its ability to assist growth, was first evidenced in numerical simulations by Ormel et al. (2007). Based on the numerical model of evolution of porosity during collisions of Suyama et al. (2008), Okuzumi et al. (2012) then investigated the formation of icy planetesimals from direct growth of porous aggre- gates and found that porosity can help grains overcome the radialdrift barrier. However, in their study, planetesimals have a density lower than 10 -3 kg m -3 , inconsistent with comet values. Kataoka et al. (2013b) thus introduced the static compression process to take into account restructuring due to gas or self-gravity compression. As it is difficult and impractical to treat at the same time the evolution of porosity during collisions at small scales and the disc physics at large scales, studies including porosity have mostly been carried out in a local context, i.e. at a given distance from the star (Ormel et al. 2007;Okuzumi et al. 2009;Kataoka et al. 2013b). However, Okuzumi et al. (2012) and Krijt et al. (2016) computed the radial evolution of grain mass and porosity in a 1D disc model. In this paper, in order to link both the small and large scales, we derive a model of porosity evolution to reproduce the effects of collisions and include it in global 3D simulations of dust growth and dynamics. We detail our models for growth and porosity as well as our numerical codes in Section 2. In Section 3, we present our results on the influence of porosity on growth and drift and we discuss on how porosity can be a solution to the radial-drift barrier in Section 4. We summarise our work and conclude in Section 5.

METHODS

Disc models

We simulate a disc of mass M disc = 0.01 M around a 1 M star. Initially, the dust-to-gas ratio by mass , defined as the ratio between dust and gas local volume densities ρ d /ρ g , is taken to be 0.01. We model the disc as power laws of the distance from the star R. Thus, the gas surface density scales as Σ g ∝ R -p and its temperature as T g ∝ R -q . In all simulations, we use a constant Shakura & Sunyaev (1973) α parameter of 10 -2 . The disc extends from R in to R out . We assume that dust moving outside of a distance R esc or inside R in is respectively escaped or accreted into the star.

We consider two disc models characterised by two different disc geometries: a disc around a classical T-Tauri star (CTTS) (studied by, e.g., Barrière-Fouchet et al. 2005;Laibe et al. 2008) and a "Flat" disc, a simple model with a constant surface density (used in, e.g., Paardekooper & Mellema 2004;Fouchet et al. 2007;Gonzalez et al. 2015a). Both models correspond to different outcomes regarding the radial-drift barrier for non-growing grains according to Laibe et al. (2012): dust grains are expected to survive the barrier in the CTTS disc but not in the Flat disc. We study both in this work to examine whether porosity impacts this behaviour. The different parameters of those disc are given in Table 1.

Drag regimes and drift

The differential velocity between gas and dust creates an aerodynamic friction. Grains are in different drag regimes according to how their size s compares to the mean free path of gas molecules λ. Stepinski & Valageas (1996) give the expression of the gas dynamical viscosity as a function of λ as

µ g = 1 2 ρ g λ c g , (1)
where

c g = k B T g m g (2)
is the gas sound speed, with m g the mean gas molecule mass. This yields

λ = 2µ g ρ g c g . (3)
The Chapman-Enskog theory provides a refined calculation of µ g . In the Sutherland model, which describes rigid elastic spheres with weak mutual attraction, it can be approximated as [START_REF] Chapman | The mathematical theory of non-uniform gases. an account of the kinetic theory of viscosity[END_REF], where σ mol = 2 × 10 -19 m 2 is the cross-section of the H 2 molecule, and we take m g = 2.32 m H = 3.883038752 × 10 -27 kg. When s < 9λ/4, the grain is in the Epstein drag regime (Epstein 1924) and when s > 9λ/4, it is in the Stokes regime [START_REF] Whipple | From Plasma to Planet[END_REF]. In those regimes, the grain's Stokes number can be expressed as

µ g = 5 √ π m g c g 64 σ mol (4)
St =                  ρ s φ s ρ g c g Ω K , (Epstein) 2 ρ s φ s 2 9 µ g Ω K , (Stokes)
(5)

where ρ s is the bulk density of the grain and Ω K is the Keplerian angular velocity. φ is the grain filling factor and will be discussed in Section 2.4. Note that some authors use the gas sound speed while others use the mean thermal velocity of gas molecules in the calculation of the drag force in the Epstein regime. We use the former. The resulting expressions for St differ by a factor √ π/8, of order unity. In the Stokes regime, even though the Stokes number would depend on the Reynolds number Re [START_REF] Whipple | From Plasma to Planet[END_REF], we deliberately do not use the Re dependency of the Stokes number for large Re for the sake of simplicity as it would require iterations to calculate it (Re depends on the differential velocity between gas and dust, which in turn depends on the stopping time, which is a function of Re). This amounts to limiting ourselves to the linear Stokes drag. Similarly, Okuzumi et al. (2012) neglected the high-Re domain of the Stokes regime (also known as the Newton drag regime) for simplicity. They showed that this somewhat accelerates dust growth but concluded that it has little effect of the ability of porosity to help grains to survive the radial-drift barrier. Additionally, the transition to the Newton regime would only occur in the inner disc regions for bodies larger than ∼ 100 m, for which selfgravity would start to be important and assist their growth. However, this is out of the scope of this paper (see Section 2.5.2).

This friction force makes the dust settle down to the mid-plane and drift towards the star. Nakagawa et al. (1986) gave the dust radial velocity as

v d, R = - St (1 + ) 2 + St 2 η v K , (6)
where v K is the Keplerian velocity and η the sub-Keplerian parameter defined as

η = - H g R 2 ∂ ln P g ∂ ln R , (7)
with H g the gas scale height and P g the gas pressure. From equation (6), we can find again that both small and large grains, with small and large St values, drift slowly. We can also infer that dust can be slowed down by increasing the dust-to-gas ratio .

Growth model

In this study, we only consider icy dust with ρ s = 917 kg m -3 . Indeed, water ice can be found in large quantities exterior to the snow line [START_REF] Van Dishoeck | Protostars and Planets VI[END_REF]. Furthermore, in those regions, water can condense on grains surface, change the sticking properties and promote growth over bouncing or fragmentation (Gundlach et al. 2011;Gundlach & Blum 2015). Silicate grains have a fragmentation threshold ∼ 1 m s -1 while icy dust fragments for a relative velocity larger than several 10 m s -1 (Blum & Wurm 2008;[START_REF] Yamamoto | [END_REF]. While recent laboratory experiments (Gundlach et al. 2018;Musiolik & Wurm 2019) found lower surface energies for water ice than previously used at low temperatures, thus suggesting that ice is not more resistant than silicates, they disagree with the tensile strengths measured in numerical simulations by Tatsuuma et al. (2019). Further investigation is thus needed before the issue can be settled. Thus, in this work we focus on pure growth and neglect bouncing and fragmentation, similarly to Krijt et al. (2016). This also allows to better understand the influence of porosity on these physical processes separately. The actual value of the fragmentation threshold of water ice therefore does not affect the results presented here. Bouncing and fragmentation will be taken into account in a followup paper.

We use a grain growth model with a locally mono-disperse mass distribution, in which collisions are considered to occur between identical grains. Details on the model implementation and a discussion on its assumptions can be found in Laibe et al. (2008). As fragmentation is not taken into account in this study, every collision leads to sticking. Stepinski & Valageas (1997) gave the expression of the variation of a grain's mass with time as it doubles in a mean collision time τ coll , i.e.

dm dt ≈ m τ coll = 4πs 2 v rel ρ d = 4πs 2 v rel ρ g . (8)
Grains collide with a relative velocity v rel transmitted from gas turbulent motion to the dust by aerodynamic drag. We use the Stepinski & Valageas (1997) model, in which

v rel = 2 3/2 Ro α √ Sc -1 Sc c g , ( 9
)
where Ro is the Rossby number for turbulent motions, taken equal to 3, and Sc the Schmidt number of the grains, whose expression is

Sc = (1 + St) 1 + ∆v 2 v 2 t , (10)
where v t = √ 2 1/2 Ro α c g denotes the turbulent velocity. ∆v = v dv g is the differential velocity between dust and gas. In our simulations, ∆v appeared to be negligible compared to v t (see also Laibe et al. 2008) and v rel can be approximated as

v rel 2 3/2 Ro α √ St 1 + St c g . ( 11
)
The reader is referred to Laibe et al. (2008) for a discussion of the difference between equation ( 10) and the expression of Youdin & Lithwick (2007), and to Laibe (2014) for a more general discussion on various models for relative velocities.

The collisional kinetic energy E kin for identical grains is then given by

E kin = 1 2 m * v 2 rel = 1 4 m v 2 rel , (12)
where the reduced mass m * = m/2 for identical grains.

Model for porosity evolution during growth

We consider grains as aggregates of elementary monomers. Those monomers are compact spheres, with a density ρ s and a radius a 0 . The volume filling factor φ is then used to characterise how porous the grains are. It is defined as the ratio between the volume occupied by matter V mat to the total volume V of a grain:

φ = V mat V = ρ ρ s ( 13
)
where ρ is the density of the grain. Compact grains have φ = 1 while φ ∼ 0 corresponds to very fluffy aggregates.

Real porous aggregates can have arbitrary shapes and mass distributions. Parameters such as radius, cross-section or volume are therefore ill-defined. The equations presented in Sections 2.2 and 2.3 were initially derived for compact spherical grains. In this work, in order to be able to use them, as well as equation ( 13), we assume that a porous aggregate can be represented by an equivalent spherically-symmetric and non-fractal collection of monomers in order to describe its spatial, size and porosity evolution.

Collisional evolution

Small grains, with low velocities [START_REF] Weidenschilling | Protostars and Planets III[END_REF] meet slowly and stick, it is the "hit-and-stick" regime. Such collisions make the grains more porous, or decrease their filling factor, as their total volume grows faster than that occupied by matter (see Okuzumi et al. 2012). When grains collide with a kinetic energy larger than the rolling energy E roll (see equation A4), i.e. the energy needed for one monomer to roll over 90 • on the surface of another monomer (Dominik & Tielens 1997), the formed aggregate dissipates the extra energy by compressing. Thus, the grain still becomes more porous, but less than in the hit-and-stick regime (Okuzumi et al. 2012). Thanks to numerical simulations of collisions, Suyama et al. (2008) have expressed the filling factor of the grain after a collision φ f in those two regimes as a function of parameters before the collision such as the kinetic energy E kin or the filling factor φ i of colliding grains. However, this formula is recursive and cannot be used in our simulations where a time step corresponds to any number (not necessarily an integer) of collisions. Thus, to describe the evolution of the filling factor during collisional growth φ col , we needed to modify the Suyama et al. (2008) model to make it continuous and non-recursive. This amounts to obtaining an expression for φ f only as a function of grain mass m and disc quantities. The details of our calculations and approximations are presented in Appendix A. Our resulting expressions for φ col are:

(i) hit-and-stick

φ h&s = m m 0 -0.58 , (14)
where m 0 = 4 3 πa 3 0 ρ s is the monomer mass. (ii) collisional compression: four cases need to be considered depending on the drag regime, Epstein or Stokes, and whether St <

1 or St > 1 φ Ep-St<1 = 2 1/5 -β Ep -3/8 β -5/8 Ep 2 1/8 × 3 10 (3/5) 5 2 3/2 Ro α m 0 c g ρ s a 0 Ω K 8 ρ g b E roll 3/8 m m 0 -1/8 , (15)
φ St-St<1 = 2 1/5 -β St -1/3 β -2/3 St × 3 10 (3/5) 5 2 3/2 Ro α m 0 c 2 g ρ s a 2 0 Ω K 36 µ g b E roll 1/3 , ( 16
)
φ Ep-St>1 = φ Ep-St<1 (M 4 ) m M 4 -1/5 , (17)
φ St-St>1 = φ St-St<1 (M 5 ) m M 5 -1/5 , (18)
where M 4 and M 5 are the masses for which St = 1, respectively in the Epstein and Stokes regimes. They are given in equations ( A16) and (A17).

Static compression

Kataoka et al. ( 2013a) have shown that very fluffy grains can be statically compressed by an applied compressive strength P. The filling factor is then linked to P as

φ = a 3 0 E roll P 1/3 . ( 19
)
Thus, the gas drag force F drag applies a compressive strength P drag = F drag /(π s 2 ) that can compact the grain up to a value φ gas depending on the drag regime:

(iii) gas drag compression

φ gas =                      m 0 ∆v π E roll ρ g c g ρ s 1/3 , (Epstein) 6 a 2 0 ∆v µ g E roll 3/8 m m 0 -1/8 . (Stokes) (20)
Finally, massive enough grains can be compressed by their selfgravity up to a filling factor φ grav given by Kataoka et al. (2013b) as:

(iv) self-gravity compression

φ grav = G m 2 0 π a 0 E roll 3/5 m m 0 2/5 . ( 21
)
The filling factor of dust naturally evolves as φ col during collisions until the grain is fluffy enough to be compressed by the gas drag or its self-gravity and then its filling factor is equal to φ gas or φ grav respectively. Therefore, the maximum of φ col , φ gas and φ grav provides the smallest possible value of the filling factor. The different regimes encountered by growing grains at fixed locations are summarised in Figure B1.

Numerical simulations

In order to model the evolution of porous dust grains in protoplanetary discs, we use two codes with different properties.

The PACED code

We have developed the PACED1 code, a 1D C++ code to model the behaviour of a single grain at a time in a stationary gas disc described with power laws with a constant dust-to-gas ratio of 10 -2 . This code contains radial drift with equation ( 6), the growth model described in Section 2.3 and porosity evolution with the model presented in Section 2.4.

3D SPH code

We also use the 3D, two-phase (gas+dust), SPH (Smoothed Particle Hydrodynamics) code described in Barrière-Fouchet et al. (2005). The code computes the dynamics of both phases, which interact via aerodynamic drag in the Epstein regime, including the backreaction of dust on gas. (In contrast to PACED, this code does not use equation ( 6) but calculates the forces on each SPH particle and their ensuing motion.) It reproduces the properties of turbulence (Arena & Gonzalez 2013). Grain growth was implemented in Laibe et al. (2008) and fragmentation in Gonzalez et al. (2015a). Selfgravity is not included. Detailed discussions of the code properties can be found in Fouchet et al. (2007) and Gonzalez et al. (2017a). In this work, we have added the implementation of the Stokes drag regime to consider the evolution of the largest solids. We have also included our model of porosity evolution.

We start the simulations presented here with 200,000 SPH gas particles arranged to reproduce the power laws presented in Section 2.1 and let the disc relax for 24 orbits at R = 100 AU for the CTTS disc and R = 40 AU in the Flat disc. Then, the same number of SPH dust particles is injected at the same locations and velocities as the gas particles, to recreate an initial dust-to-gas ratio = 10 -2 . Both sets of particles are then evolved together. Previous work with the same code (e.g. Barrière-Fouchet et al. 2005;Fouchet et al. 2007;Pignatale et al. 2019) has shown that convergence is reached at lower particle numbers. We also verify that the resolution criterion from Laibe & Price (2012) is met. Each SPH dust particle represents a collection of identical physical dust grains (locally mono-disperse assumption, see Section 2.3) and carries intrinsic dust properties such as mass, filling factor and thus, size, that are described according to the models presented in Sections 2.3 and 2.4. In particular, we chose a uniform initial size s 0 = 10 µm for all grains to shorten the computation times. Indeed, we have previously shown that very small grains grow fast and quickly forget their initial size (Laibe et al. 2008). Their corresponding initial filling factor and mass, which depend on their location, are precalculated with the PACED code by evolving compact monomers of size a 0 without radial drift. The simulations are stopped when the largest grains reach kilometre sizes, for which self-gravity becomes important, because both our growth model and our SPH code do not take it into account. Obtaining estimates of the fraction of dust that would turn into planetesimals or of properties of the resulting planetesimal population is therefore not feasible in the framework of our study.

Comparison

As summarised in Table 2 and contrary to the PACED code, our 3D SPH code models the evolution of the gas disc and takes into account the collective effects. Indeed, dust drift can change the local dust-to-gas ratio which has an impact on dust growth and dynamics, while it is taken constant in the PACED code. However, the 3D SPH code requires long calculation times and is not practical to study a large variety of cases. Thus, we use the PACED code to test our model and estimate the dust behaviour on larger time scales. Moreover, in the PACED code, disc quantities such as gas density or differential velocity between gas and dust are obtained from analytical prescriptions while in the 3D SPH code, they are computed self-consistently (Barrière- Fouchet et al. 2005).

RESULTS

In Sections 3.1 and 3.2, we present how porosity evolves during dust growth and how it impacts growth and dynamics in the CTTS disc with monomer size a 0 = 0.1 µm. In Section 3.3, we show the influence of the monomer size a 0 on porosity and thus, on growth in the Flat disc.

Evolution of the filling factor

Figure 1 shows the evolution of the filling factor of grains of different initial positions as a function of their mass as they grow and drift radially, computed with the PACED code. The curves have similar shapes to those found by Kataoka et al. (2013b) for grains at fixed distances from the star (see also our model at fixed distances in Fig. B1 for comparison): φ first decreases due to collisions, then increases again with gas drag compression before selfgravity is strong enough to compact the grains further. However, inwards drifting grains experience different disc conditions, and in particular larger gas densities, than static ones and are all the more compacted by gas drag as they are close to the star. Their minimum filling factor, ∼ 10 4 , is somewhat larger and it is reached for smaller masses.

The time evolution of the grain filling factor computed with the 3D SPH code is displayed in Fig. 2. Initially, 10-µm sized grains have a filling factor of ∼ 10 -1 -10 -2 , according to their initial location in the disc. As dust grows, representative dots move to the right. At first, grains are in the collisional compression regime. This regime tends to make grains more and more porous until they are fluffy enough to be compressed by gas drag. Note that the gas drag compression regime has a different shape than that computed with the PACED code (Fig. 1), illustrating the differences between both codes listed in Section 2.5.3. In the top left panel, some dots move upwards with a quasi constant mass. Those dots correspond to grains that are quickly compacted by gas drag as they are accreted into the central star. The minimum value of φ, also ∼ 10 4 , is consistent with the numerical simulations of collisions of Kempf et al. (1999).

Figure 3 shows the distribution of the filling factor through the disc. A vertical sorting of the filling factor occurs: more compact grains are located along the disc surface while grains are fluffier and fluffier as they are closer to the mid-plane. Indeed, grains grow more rapidly in denser regions (equation 8) and thus reach lower filling factors in the collisional regime. The same effect is seen in the radial direction. However, grains in the disc innermost regions are compacted because it is there that gas drag compression is mostly effective, where the gas is dense enough.

Influence on grain size and dynamics

Figure 4 compares the radial evolution of the size of compact and porous grains, computed with the PACED code. Compact and porous grains both have a qualitatively similar behaviour: they first grow at a quasi-constant distance from the star until they approach an optimal size s St=1 , for which St = 1 and radial drift is the fastest. They then drift rapidly inwards while continuing to grow before slowing down close to the star, where they finally keep growing while drifting very little. These three stages were described and ex- plained by Laibe et al. (2008). During the second stage, compact grains gain 1 to 3 orders of magnitude in size before reaching the disc inner regions, depending on their initial location, while porous grains can gain up to 4. Indeed, since s St=1 is larger for porous grains (see Fig. 9), they have a larger cross section and grow more rapidly (equation 8): the slope in Fig. 4 is steepest. These results are in agreement with those found by Okuzumi et al. (2012).

The time evolution of the radial grain size distribution of compact and porous grains computed with the 3D SPH code is plotted on Fig. 5. Here as well, porous grains experience a quicker growth than compact grains. In 8090 yr, in the inner disc, the largest porous grains have planetesimal sizes and are in the Stokes drag regime while compact grains hardly reach one meter and stay in the Epstein regime. Furthermore, growth is very slow for compact dust beyond 100 au and grains cannot grow beyond 100 µm. On the contrary, porous dust can reach centimetre sizes up to 300 au and metre sizes interior to 100 au. The right column of Fig. 5 shows that small grains are only slightly porous. As they grow, they become fluffier and fluffier and are compressed in the inner regions of the disc as seen in Fig. 3 as well.

Fig. 6 shows the time evolution of the dust spatial distribution of compact and porous grains computed with the 3D SPH code. In both cases, dust is vertically size-sorted. Larger solids can be found close to the mid-plane while small grains are distributed over a larger scale height. In the porous case, the mid-plane is filled with ∼ 100 m solids from the inner edge out to ∼ 50 au. Two differences on dynamics can be also spotted out. After 8090 yr, the disc made of compact dust is approximately 20 au less radially extended than with porous grains. Moreover, this disc is also more settled compared to the one with fluffy dust. Indeed, porous grains can stay coupled with the gas phase (St 1) at larger sizes than compact grains, they settle down and drift inward slightly less rapidly.

Influence of monomer size

Figure 7 compares the time evolution of the filling factor during growth for grains made of monomers of different sizes a 0 = 0.1 and 1 µm in the Flat disc. The behaviours are similar to the CTTS disc (Fig. 2). Initially, grains are in the collisional expansion regime then are compacted because of the gas drag. Grains reach similar masses in both cases, large enough to be compressed by their selfgravity as indicated by the red arrow. All these regimes depend on the monomer size and the filling factor increases when a 0 increases, see equations ( 14)-( 21). Thus, 10-µm sized grains composed of 1-µm monomers are more massive than those made of 0.1-µm monomers. One can observe that gas drag compression is almost as efficient in both cases, i.e. it compresses grains to similar filling factors. However, since grains made of large monomers are already more compact, their change in filling factor is smaller during that phase. In the same way, self-gravity starts to compact the dust at a larger mass.

The radial size distribution is compared in Fig. 8. At first, the difference in filling factor does not affect much the radial extension of the disc, i.e. the radial drift is not influenced by the monomer size. The main impact of the porosity change is on growth. Since grains composed of large monomers are less porous, they grow less efficiently and produce less large grains. Dust made of 0.1µm monomers can produce planetesimals in the innermost 100 au of the disc while it is only the case over a few au with of 1-µm monomers. 2012) studied the possibility of forming planetesimals from direct growth of porous dust. This idea was then completed by Kataoka et al. (2013b) who added the compression due to gas drag and self-gravity of massive solids. They have thus determined the evolution of the filling factor during growth but without drift and collective effects in a Minimum Mass Solar Nebula (MMSN) disc. With our continuous model, the values we obtain are generally in a good agreement with those from Kataoka et al. (2013b). The minimum value reached by the filling factor, ∼ 10 -4 , is however slightly larger than that of Kataoka et al. (2013b) because our grains drift to regions where the gas density is larger, and hence gas drag compression stronger.

The Solar System still contains porous bodies today. Evidence has been found for low densities in asteroids, such as 253 Mathilde (Housen et al. 1999;Veverka et al. 1999), or meteorites (e.g. Brown et al. 2002). This may even be a more widespread property, with the recent suggestion that the first known interstellar object to visit our system, 1I/'Oumuamua, is a fractal dust aggregate with an ultra-low density of 10 -2 kg m -3 (Flekkøy et al. 2019). Spatial exploration and measurements on comets have shown that these objects made of ice and silicates are also porous. Their average density falls between 400 and 600 kg m -3 (Blum et al. 2006;A'Hearn 2011). As an example, Sierks et al. (2015) determined that the density of comet 67P/Churyumov-Gerasimenko is of the order of 470 kg m -3 . The Rosetta probe orbiting comet 67P and its lander Philae have provided a wealth of data on this body. In particular, they have shown that its surface is covered with solids whose size varies from tens of micrometers to tens of meters (Blum et al. 2017). Among them, millimiter-sized grains with a density < 1 kg m -3 and a filling factor φ ∼ 10 -3 were found with the GIADA instrument (Fulle et al. 2015). As can be seen on Fig. 5, our results show that grains with similar sizes and filling factors are present over large regions in the disc. While such low filling factors can seem surprising at first, they are compatible with the measurements made on comet 67P.

The Stokes regime, a key point

One of the main theoretical problems for dust evolution is the radial-drift barrier, i.e. the rapid inwards drift of grains with St = 1 leading to dust accretion onto the star. A solution to this problem is either to stop the grains from drifting or for them to grow from St 1 to St 1 (for which the drift is slow) in a time shorter than the drift timescale.

Moreover, Laibe et al. (2012) have shown that if dust reaches St > 1 in the Stokes drag regime, it remains in the disc if q ≤ 2/3. Discs satisfying this condition represent 90% of discs observed by Andrews & Williams (2005, 2007). The Stokes regime has two advantages: the Stokes number increases as grains get closer to the star (while it decreases in the Epstein regime) and varies as s 2 (while as s in the Epstein regime), see equation ( 5) 2 . Thus, in the Stokes drag regime, dust can reach large Stokes numbers more easily. Physically, it means that the grain is more efficiently slowed down by the stronger gas drag in the inner parts of the disc. Consequently, reaching the Stokes drag regime with St ∼ 1 is a key point 2 For a power-law disc, St ∝ s R p in the Epstein regime and St ∝ s 2 R (q-3)/2 in the Stokes regime. for grains to remain in most discs and thus survive the radial-drift barrier.

In order to understand the influence of porosity on the Stokes number and the drag regime, one needs to compare two characteristic sizes:

• the optimal size, s St=1 , for which the grain Stokes number reaches unity, i.e. the size corresponding to the fastest drift;

• 9λ/4, the transition size between the Epstein and Stokes drag regimes.

In a power-law disc, the optimal size s St=1 ∝ φ -1 R -p in the Epstein regime and s St=1 ∝ φ -1/2 R (3-q)/4 in the Stokes regime. 9λ/4 does not depend on the filling factor and varies as R p+(3-q)/2 . If s St=1 < 9λ/4 (resp. s St=1 > 9λ/4), grains reach the maximal drift velocity in the Epstein (resp. Stokes) regime.

Figure 9 shows the optimal size for different filling factors and the transition between both drag regimes. For compact grains, the transition between the Epstein and Stokes regimes occurs necessarily with St > 1. It means that compact grains experience a maximal inwards drift (and accretion onto the star) without the chance to transition to the Stokes regime. On the contrary, for grains with a filling factor lower than 0.5, there exists a region in the inner disc where porous grains are in the Stokes drag regime with St ≤ 1, which is more and more extended as grains become fluffier and fluffier. Thus, the Stokes number of those grains can increase rapidly during their growth, allowing them to decouple from the gas and survive the radial-drift barrier.

Porosity and growth

As discussed previously, porosity can have an impact on dust spatial evolution and growth. According to Figs. 4 and5, compact and porous grains have a qualitatively similar growth behaviour in three stages. The optimal size for compact grains is about 0.1 mm. Since the optimal size depends on the grain filling factor, porous grains can reach larger sizes before starting to drift rapidly. For instance, porous grains with φ = 10 -3 , order of magnitude of the average filling factor over the disc (Fig. 5), have an optimal size 10 3 times larger than for compact grains. Thus, they can grow up to ∼ 10 cm before drifting significantly. This average ratio of 10 3 for porous grains sizes can be seen in Figs. 4 and5. Generally, porous dust is characterised by a more efficient growth during the three steps described previously as their collision cross-section is larger than for compact grains. Consequently, porosity helps grains to reach larger sizes and transition to the Stokes drag regime in the inner disc, as seen in Figs. 4 and5. We discussed in Section 4.2 on how the Stokes drag regime is important for the dust survival in the disc.

Porosity also allows to form 100 m -1 km planetesimals in the mid-plane in the innermost 100 au, as shown by Figs. 4 and6, in less than 10 4 yr. Even if our disc model is different from the one Krijt et al. (2016) use, we find drift and growth timescales for both compact and porous grains similar to theirs. Those planetesimals are thought to be building blocks to form giant planet cores and our results can be used as a starting point for simulations of giant planet formation (Kobayashi et al. 2016;Kobayashi & Tanaka 2018).

On drift and the importance of collective effects

We can spot a different behaviour for compact dust in the very inner disc. With the PACED code, compact dust is accreted onto the star without growing larger than the decimetre (Fig. 4), while we find in 3D SPH simulations that compact solids can remain in the very inner disc where they can reach metre sizes (Fig. 5). This phenomenon can be explained by the different calculation (and influence) of the local dust-to-gas ratio in our two codes. Indeed, Nakagawa et al. (1986) have shown (via equation ( 6)) that dust drift is slowed down when the local dust-to-gas ratio increases thanks to dust collective effects. Moreover, the PACED code uses a constant of 10 -2 while in the 3D SPH code, this quantity is directly calculated using the local dust and gas densities and so may vary. Figure 10 reports the value of the dust-to-gas ratio at several times of the simulation with our 3D SPH code. Initially, = 10 -2 in the whole disc as in the PACED code but as grains grow, they drift inwards and the dust-to-gas ratio increases in the innermost few au where its vertically integrated value can reach 10 -1 . In the midplane, it can even exceed unity thanks to vertical settling. Thus, the drift is slowed down and growth accelerated as the dust phase is denser. Those two effects help compact grains to reach meter sizes at the edge of the disc. However, they are not strong enough to allow them to reach the Stokes regime.

Such an increase in the dust-to-gas ratio in the inner disc was not seen for compact grains by Birnstiel et al. (2010), Okuzumi et al. (2012), or Krijt et al. (2016), whose framework and assumptions are very different from ours. These previous studies are all one-dimensional and compute the vertically-averaged dust evolution in the radial direction while keeping the gas distribution fixed (amounting to neglecting the back-reaction of dust on gas) and assuming vertical hydrostatic equilibrium for the gas and a prescription for the dust vertical scale-height as a function of turbulence. Okuzumi et al. (2012) further assess the impact of back-reaction on the velocities of gas and dust and on collision velocities and find it to be negligible. However, it is not a full implementation of backreaction as the background gas structure is still kept constant. For compact dust, all these studies found that grains grow until they reach St = 1 then drift inwards, producing a drift-limited size distribution and a flat dust-to-gas ratio profile. On the contrary, our simulations are three-dimensional and compute the self-consistent evolution of both gas and dust, including back-reaction of dust on gas. It has been shown that back-reaction, assisted by an increase of dust density close to the mid-plane due to vertical settling, is able to slow down the dust radial drift enough to allow grains to grow to sizes for which St > 1, slow down their drift further, and pile up (e.g. Laibe et al. 2008;Gonzalez et al. 2015bGonzalez et al. , 2017a)), something models neglecting back-reaction are unable to capture. The importance of back-reaction has further been shown by several authors (e.g. Kanagawa et al. 2017;Dipierro & Laibe 2017;Dipierro et al. 2018). The reader is referred to Gonzalez et al. (2017a) for a comparison of dust-to-gas ratio profiles with and without back-reaction in the case when dust fragmentation is included.

For porous dust, Okuzumi et al. (2012) and Krijt et al. (2016) found that grains grow rapidly, overcome the radial-drift barrier, then pile up, leading to an increase of the dust-to-gas ratio in the inner disc. This is exactly what we find in this study, see Section 3.2 and Fig. 10. However, where Okuzumi et al. (2012) saw an increase of their one-dimensional of a factor of several, our azimuthally-averaged, vertically-integraged is enhanced by one order of magnitude. Porous grains, similarly to compact grains, experience collective effects, which play an important role in setting their spatial distribution. Additionally, even though porosity can slighty slow down vertical settling and radial drift, its main effect remains a strong acceleration of growth (Section 3.2), which causes the transition to the Stokes drag regime. As a result, porous grains decouple from the gas and pile-up sooner in their evolution, therefore at larger distances from the star, than compact grains. This is reflected in the locations of the dust-to-gas ratio maxima in Fig. 10. The combination of collective effects and of the Stokes regime also results in a larger dust-to-gas ratio than when dust back-reaction is not included.

CONCLUSION

The growth from sub-µm monomers to planetesimals is hampered by several barriers such as the radial-drift barrier. The dust needs to decouple from the gas in order to remain in the disc. In this work, we investigate how porosity can act on both grain growth and drift to overcome this barrier. It amounts to considering simultaneously, and thus link, the small (collisions) and large (disc) scales. To do so, we have developed a model of evolution of the filling factor depending on grain characteristics (bulk density, rolling energy, mass) and disc quantities (gas density, sound speed or temperature). Our model can be used for codes in which the disc evolution time step is different from that of individual collisions. We demonstrate that grains can remain in the disc in most cases if they grow enough to transition to the Stokes drag regime. Indeed, in this regime, dust can reach large Stokes numbers and decouple from the gas more easily. We have shown that compact grains do not grow quickly enough and stay in the Epstein regime. However, collective effects can slow them down enough for them to remain in the very inner regions of the CTTS disc. We find that the growth is accelerated for porous grains, allowing them to transition to the Stokes regime close to the star and survive the radial-drift barrier, both in the CTTS and Flat discs. Furthermore, our study has shown that porous millimetre grains have an average filling factor of about 10 -3 , in good agreement with measurements made on comet 67P/Churyumov-Gerasimenko (Fulle et al. 2015). Finally, we find that small planetesimals can be formed by direct coagulation of porous dust in the innermost 100 au. This result provides a link with the formation of giant planet cores from planetesimals.

       (3/5) 5 (E kin -3 b E roll ) N 5 tot b E roll V 10/3 0 + 2V 5/6 i -4        3/10 , for E kin ≥ 3 b E roll (internal restructuring regime) (A1)
where m i is the mass of colliding grains and V i their volume. N tot is the total number of monomers involved in the collision, V 0 is the volume of a monomer, defined as a compact sphere of radius a 0 and b is a numerical factor taken as 0.15 (Okuzumi et al. 2012).

The rolling energy is given by Dominik & Tielens (1997) as

E roll = 6π 2 γ a 0 ξ crit , ( A2
)
where γ is the surface energy of the material. They find that ξ crit , the critical rolling distance of one monomer on another for energy dissipation, is of the same order as the critical distance for monomer separation, given by [START_REF] Chokshi | [END_REF] as

δ c = 27π 2 γ 2 a 0 2 E 2 1/3 , ( A3
)
where E is the Young modulus of the material. This leads to

E roll = 2916 π 8 γ 5 a 4 0 E 2 1/3 , ( A4
)
where we adopt for ice γ = 7.3 × 10 -2 J m -2 and E = 9.4 GPa [START_REF] Yamamoto | [END_REF].

In this Appendix, we explain how we extend the Suyama et al. ( 2008) equation (A1) to make it continuous in Section A1 for the hit and stick regime (i.e. E kin ≤ 3 b E roll ) and in Section A2 for the internal restructuring regime (i.e. E kin ≥ 3 b E roll ). It is equivalent to express φ f as a power law of the grain mass m and quantities of the disc. To do so, we chose to approximate E kin as either very small or very large compared to E roll in order to do a finite expansion. More detailed calculations leading to the equations presented in this Appendix can be found in Garcia (2018, in French).

A1 Hit-and-stick regime

If we consider E kin E roll , equation (A1) becomes

φ f = 2 2.99 φ i . (A5)
Thus, we see that the filling factor in the hit-and-stick regime evolves as a geometrical progression with a common ratio 2/2.99. We can express the filling factor φ f of a grain with a mass m as a function of the monomer filling factor φ 0 and n the number of collisions to form that grain

φ f = φ 0 2 2.99 n .
(A6)

After one collision, the grain mass doubles. So after n collisions from a monomer, m = 2 n m 0 . Consequently, the filling factor φ f in the hit-and-stick regime (thereafter renamed φ h&s ) can be given as a function of m φ h&s = φ 0 m m 0 ln(2/2.99)/ ln(2)

.

(A7)

Since monomers are compact, φ 0 = 1. Note that we no longer have the recursive aspect of equation (A1).

A2 Internal restructuring regime

We consider here that E kin E roll . Using N tot = 2m i /m 0 and V i /V 0 = 1/φ i (m i /m 0 ), equation (A1) becomes

φ f = 2 -1/5 φ i        1 + (3/5) 5 E kin 2 b E roll 1 φ 10/3 i m i m 0 -5/3        3/10 . (A8)
However, E kin depends on φ i trough v 2 rel and the Stokes number St. Nevertheless, v 2 rel does not vary linearly with the Stokes number. We consider St 1 (resp. St 1) in order to have E kin ∝ St (resp. St -1 ). With this approximation, φ f is then related to φ i at a certain power. The power depends on the grain drag regime (Epstein or Stokes) and if its Stokes number is smaller or larger than 1. As we want to express φ f as φ f ∝ m k with k ∈ R, we have for given disc quantities, φ f = 2 k φ i . We define β = 2 -k , and thus φ i = βφ f . The values of β are different in Epstein and Stokes regimes and are taken to fit Eq. (A1) as discussed in Sections A2.1 and A2.2. codes. Our model deviates from that of Suyama et al. (2008) in two situations: for E cin ∼ E roll and St ∼ 1. In both cases, the filling factor value we obtain is slightly higher than obtained with the discrete model but this difference is not significant for grains growth and dynamics and can be neglected.

A4 Algorithm

The algorithm allowing to compute the value of the filling factor after collisions φ col is detailed in Algorithm 1.

Figure 1 .

 1 Figure 1. Filling factor as a function of mass in the CTTS disc for drifting grains of different initial positions R 0 with the PACED code. Arrows show the different regimes encountered by the grains.

Figure 2 .

 2 Figure 2. Time evolution of the filling factor as a function of mass for grains in the CTTS disc with the 3D SPH code. Representative dots move to the right as grains grow. The filling factor decreases thanks to the collisional expansion regime (green arrow) then increases because of the compression due to the gas drag (blue arrow).

Figure 3 .

 3 Figure 3. Time evolution of the dust distribution in the CTTS disc with the 3D SPH code. The colour bar represents the filling factor.

Figure 4 .

 4 Figure 4. Grain size as a function of the distance from the star in the CTTS disc for different initial positions R 0 with the PACED code, for compact (left) and porous (right) dust. The red dashed line marks the limit between the Epstein and Stokes regimes. The symbols show the value of the size at the same time steps as in Fig.5.

Figure 5 .

 5 Figure 5. Radial grain size distribution in the CTTS disc obtained with the 3D SPH code. Porous grains (right) can reach larger sizes than compact dust (left). The colour bar represents the filling factor. The red dashed line in the top left corner separates the Epstein (below the line) and Stokes (above the line) drag regimes. Four snapshots at 100, 2490, 4880 and 8090 yr are shown, from top to bottom.

Figure 6 .

 6 Figure 6. Time evolution of the dust distribution in the CTTS disc obtained with the 3D SPH code. The colour bar represents the grain size. Porous grains (right) are less settled than compact dust (left). Four snapshots at 100, 2490, 4880 and 8090 yr are shown, from top to bottom.

Figure 7 .

 7 Figure 7. Time evolution of the filling factor as a function of mass for grains in the Flat disc obtained with the 3D SPH code. Porous grains that are made of 0.1 µm monomers (left) can become more porous than grains formed with 1 µm monomers (right). However, the filling factor has a similar behaviour in both cases: it decreases thanks to the collisional expansion regime (green arrow) then increases because of the compression due to the gas drag (blue arrow) and finally, grains reach the self-gravity compression regime (red arrow). Four snapshots at 44, 6470, 12900 and 19330 yr are shown, from top to bottom.

Figure 8 .

 8 Figure 8. Time evolution of the radial grain size distribution in the Flat disc obtained with the 3D SPH code. Aggregates made of 0.1-µm monomers (left) can reach larger sizes all over the disc than dust made of 1-µm monomers (right). The colour bar represents the filling factor. Four snapshots at 44, 6470, 12900 and 19330 yr are shown, from top to bottom.

  in the Solar SystemOkuzumi et al. (

Figure 9 .

 9 Figure 9. Optimal size, i.e. size for which St = 1, as a function of the distance from the star in the CTTS disc for different filling factor values. For a given filling factor, below the line, grains have St < 1 and above it St > 1. The red line represents the limit between the Epstein and Stokes regimes for s = 9λ/4. Above this line, grains are in the Stokes regime and St ∝ s 2 while below, they are in the Epstein regime and St ∝ s.

Figure 10 .

 10 Figure10. Azimuthally-averaged radial profiles of the dust-to-gas ratio: vertically integrated (top) and in the mid-plane (bottom) for different time steps obtained with a 3D SPH simulation of the CTTS disc. Porous grains are represented by solid lines and compact grains by dotted lines. In the PACED code, the dust-to-gas ratio has a constant value of 10 -2 .

A2. 1

 1 In the Epstein regime and St < 1In the Epstein regime with St < 1, the filling factor φ f (thereafter renamed φ Ep-St<1 ) of a grain with mass m can be expressed asφ Ep-St<1 = 2 1/5 -β Ep -3/8 β -5/8 Ep 2 1/8 × 3 10 (3/5) 5 2 3/2 Ro α m 0 c g ρ s a 0 Ω K 8 ρ g b E roll3Ep is the value of β in the Epstein regime. As k = -1/8, β Ep = 2 1/8 .

Figure A2 .

 A2 Figure A2. Comparison of our model of filling factor evolution during collisions (straight cyan lines) with the Suyama et al. (2008) discrete model (red dots) for grains evolving at fixed positions of 5 au (left) and 100 au (right) in the CTTS disc.

Figure B1 .

 B1 Figure B1. Grain filling factor as a function of mass in the CTTS disc at different fixed positions R 0 with the PACED code. Arrows show the different regimes encountered by the grains.

Table 1 .

 1 Parameters of our disc models. Values with subscript 0 are given for R 0 = 1 AU.

	Model	p	q	R in	R out	R esc	T g,0	Σ g,0
				(AU) (AU) (AU) (K)	(kg.m -2 )
	CTTS	3/2 3/4 3	300	400	197	4537
	Flat	0	1	4	120	160	619	11

Table 2 .

 2 Comparison between our PACED and 3D SPH codes.

	Characteristics	PACED 3D SPH
	dynamics	X	X
	growth	X	X
	porosity	X	X
	turbulence		X
	disc evolution		X
	collective effects		X

© 2020 The Authors

MNRAS 000, 1-13 (2020)

for Porous And Compact dust Evolution in Discs MNRAS 000, 1-13 (2020)

ACKNOWLEDGEMENTS

We thank the anonymous referee for their careful reading and assessment of our work, which helped to improve this paper. The authors acknowledge funding from ANR (Agence Nationale de la Recherche) of France under contract number ANR-16-CE31-0013 (Planet-Forming-Disks). This research was partially supported by the Programme National de Physique Stellaire and the Programme National de Planétologie of CNRS/INSU, France. AJLG and JFG thank the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) of the Université de Lyon for its financial support within the programme 'Investissements d'Avenir' (ANR-11-IDEX-0007) of the French government operated by the ANR. Simulations were run at the Common Computing Facility (CCF) of LABEX LIO. All figures were made with with the Python library matplotlib (Hunter 2007).

APPENDIX A: MODEL OF POROSITY EVOLUTION DURING COLLISIONS Suyama et al. (2008) have found that the filling factor after collision φ f is related to the filling factor φ i and some other quantities such as the kinetic energy E kin before collision by In the Stokes regime with St < 1, the filling factor φ f (thereafter renamed φ St-St<1 ) of a grain with mass m can be expressed as

where M 4 and M 5 are respectively given by equations ( A16) and (A17).

A3 Transition masses

In order to know in which regime the grain is, we compare its mass m with the limit masses in every regime with the algorithm. The names of those transition masses are reported in Table A1. The transition masses M 1 to M 5 are given by

Ep 36 8/3 8 -3 × 3 10

)

Note that those expressions depend only on the distance from the star R in the case of a power-law disc model. The different transition masses are plotted in Fig. A1 for the CTTS disc. As shown in Figs. A1 andA2, all the dust begins to grow in the hit-and-stick regime before getting compacted by collisions in the Epstein drag regime. However, two different cases can be highlighted: grains in the first tens of AU transition to the Stokes regime as they grow while further in the disc, they keep growing in the Epstein regime. Moreover, Fig. A2 exhibits a comparison between the Suyama et al. (2008) discrete model and our continuous model. Both models mostly give similar results, allowing the simpler implementation of the continuous model in hydrodynamical

APPENDIX B: POROSITY EVOLUTION OF STATIC GRAINS IN THE FULL MODEL

Figure B1 shows the evolution of the filling factor of static grains as they grow at different distances from the star as a function of their mass, computed with the PACED code. The different regimes in the collisional evolution (hit-and-stick and collisional compression) and static compression (gas drag and self-gravity) phases are identified. Figure B1 is to be compared with Fig. 1, showing the porosity evolution of radially drifting grains (see Section 3.1). This paper has been typeset from a T E X/L A T E X file prepared by the author.