Émilie Chouzenoux

Víctor Elvira

GRAPHEM: EM ALGORITHM FOR BLIND KALMAN FILTERING UNDER GRAPHICAL SPARSITY CONSTRAINTS

Keywords: Kalman filtering, state-space model, graphical inference, sparsity, proximal methods, EM algorithm

Modeling and inference with multivariate sequences is central in a number of signal processing applications such as acoustics, social network analysis, biomedical, and finance, to name a few. The linear-Gaussian state-space model is a common way to describe a time series through the evolution of a hidden state, with the advantage of presenting a simple inference procedure due to the celebrated Kalman filter. A fundamental question when analyzing a multivariate sequence is the search for relationships between its entries (or the entries of the modeled hidden state), especially when the inherent structure is a non-fully connected graph. In such context, graphical modeling combined with parsimony constraints allows to limit the proliferation of parameters and enables a compact data representation which is easier to interpret by the experts. In this work, we propose a novel expectation-maximization algorithm for estimating the linear matrix operator in the state equation of a linear-Gaussian state-space model. Lasso regularization is included in the M-step, that we solve using a proximal splitting Douglas-Rachford algorithm. Numerical experiments illustrate the benefits of the proposed model and inference technique, named GraphEM, over competitors relying on Granger causality.

INTRODUCTION

State-space models allow for the statistical description of systems in countless applications in science and engineering. They model a hidden state that evolves over (discrete) time, which is only partially observed through transformed and noisy partial observations. In its simple form, the state-space models have a Markovian dependency in the hidden-state, and each observation depends on the current state, being conditionally independent on the rest of states.

V.E. and É.C. acknowledge support from the Agence Nationale de la Recherche of France under PISCES (ANR-17-CE40-0031-01) and MAJIC (ANR-17-CE40-0004-01) projects.

In this framework, a common statistical task is the estimation of the hidden state conditioned to the available observations. When the estimation is done in a probabilistic manner, this is commonly known as Bayesian filtering. Similarly, when the estimation of each state is done using the whole sequence of observations, the task is called Bayesian smoothing. The linear-Gaussian state-space model has been widely used in a plethora of problems and admits an exact computation of both filtering and smoothing distribution through the Kalman filter and the Rauch-Tung-Striebel (RTS) smoother [START_REF] Sarkka | Bayesian Filtering and Smoothing[END_REF]Chapter 8]. For non-linear state-space models, there exist many other methods such as the extended KF (EKF), the unscented KF (UKF) [START_REF] Julier | Unscented filtering and nonlinear estimation[END_REF], or the particle filters (PF) [START_REF] Djurić | Density assisted particle filters for state and parameter estimation[END_REF]. However, all these algorithms need the model parameters to be known. In the Kalman filter literature, there exist methods for the estimation of the static parameters, based on the expectation-maximization (EM) algorithm or optimization-based methods [START_REF] Sarkka | Bayesian Filtering and Smoothing[END_REF]Chapter 12].

In this paper, we propose a method called GraphEM for estimating the linear matrix of the state equation of a linear-Gaussian state-space model. While this parameter is arguably the hardest to be estimated, it provides very valuable information about the hidden process, not only for inference purposes but also for understanding the uncovered relations among the state dimensions, in the line of graphical modeling methods for time series [START_REF] Eichler | Graphical modelling of multivariate time series[END_REF][START_REF] Bach | Learning graphical models for stationary time series[END_REF][START_REF] Barber | Graphical models for time-series[END_REF]. Such representation of multivariate sequences interactions has applications in several domains such as biology [START_REF] Pirayre | BRANE Clust: Cluster-assisted gene regulatory network inference refinement[END_REF][START_REF] Luengo | Hierarchical algorithms for causality retrieval in atrial fibrillation intracavitary electrograms[END_REF], social network analysis [START_REF] Ravazzi | Learning influence structure in sparse social networks[END_REF] and neurosciences [START_REF] Richiardi | Machine learning with brain graphs[END_REF].

Without any loss of generality, the GraphEM method is developed with the following setup in mind. The multidimensional state contains several uni-dimensional time series, each of them corresponding to a particular node. The matrix in the state equation encodes the linear dependency between the state of a given node and the previous states in the neighbor nodes. GraphEM aims at estimating the directed graph, i.e., the connections between nodes and the associated weights, that represents the (causal) dependencies between the states. GraphEM imposes a sparsity constraint in the linear matrix which promotes that each dimension of the state depends of only small subset of states in the previous time step. The method relies on the EM framework, adapted for maximum a posteriori (MAP) estimation thanks to the introduction of an ℓ 1 penalization function. The E-step amounts to apply the Kalman filtering and RTS smoothing recursion. The M-step aims at solving a Lasso-like problem, that we address through a Douglas-Rachford proximal splitting algorithm, benefiting from great practical performance and sounded convergence guarantees [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF]. Experimental results on several synthetic datasets allow to assess the fast practical convergence rate of our method, and its ability to estimate efficiently the causality graph, in particular when compared to state-of-the-art techniques relying on Granger causality.

The paper is structured as follows. In Section 2 we briefly describe the model and the filtering/smoothing algorithms, and introduce the EM framework for parameters estimation in state-space models. The novel GraphEM algorithm is presented in Section 3. The paper concludes with a numerical example in Section 4 and some remarks in Section 5.

BACKGROUND

Linear state-space model

Let us consider the following Markovian state-space model,

x k = Ax k-1 + q k , y k = Hx k + r k , (1)
for k = 1, . . . , K, where x k ∈ R Nx is the hidden state at time k, y k ∈ R Ny is the associated observation, A = R Nx×Nx , H = R Ny×Nx , {q k } K k=1 ∼ N (0, Q) is the i.i.d. state noise process and {r k } K k=1 ∼ N (0, R) is the i.i.d. observation noise process. The state process is initialized as x 0 ∼ N (x 0 ; x0 , P 0 ) with known x0 and P 0 .

Kalman filtering and smoothing

The interest in the linear-Gaussian model of Eq. ( 1) is usually in the computation of the sequence of filtering distributions p(x k |y 1:k ), where we use the shorthand y 1:k = {y j } k j=1 . The Kalman filter provides an exact form of those as p(x k |y 1:k ) = N (x k |m k , P k ), for every k = 1, . . . , K [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]. The sequence of smoothing distributions (conditioned to the whole set of observations), can be also computed exactly by the RTS smoother yielding p(x k |y

1:K ) = N (x k |m s k , P s k ).
In both Kalman filter and RTS smoother, it is necessary to know the exact model, i.e., the parameters A, H, Q, and R. The explicit filter and the smoother algorithms can be found in [START_REF] Sarkka | Bayesian Filtering and Smoothing[END_REF].

EM framework for parameters estimation

In this work, we are interested in the joint estimation of the unknown A and the hidden states. Parametric estimation in state-space models has been considered in the literature through three main types of methods: EM algorithms [START_REF] Shumway | An approach to time series smoothing and forecasting using the EM algorithm[END_REF][START_REF] Thiesson | Arma time-series modeling with graphical models[END_REF], optimization based methods [START_REF] Olsson | State-space models: from the EM algorithm to a gradient approach[END_REF], and Monte Carlo methods [START_REF] Kantas | An overview of sequential Monte Carlo methods for parameter estimation in general state-space models[END_REF]. We will retain here the EM-based framework since it is well-suited for linear Gaussian models and it is flexible while benefiting from sounded convergence guarantees. EM is a method to iteratively find a maximum likelihood (ML) estimate of the parameters when the direct optimization of the posterior distribution is not feasible [START_REF] Moon | The expectation-maximization algorithm[END_REF][START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF]. It is also highly connected to majoration-minimization (MM) approaches, as it can be viewed as a special class for those [START_REF] Hunter | A tutorial on mm algorithms[END_REF]. The EM algorithm indeed alternates between a majoration step consisting in building an upper bound on the neg-loglikelihood function (E-step), and the minimization of this upper bound (M-step). In the context of state-space models with parameters θ, at each iteration i ∈ N of EM method, one majorizeslog p(y 1:K |θ) ≤ Q(θ; θ (i-1) ) [START_REF] Shumway | An approach to time series smoothing and forecasting using the EM algorithm[END_REF] with θ (i-1) the parameters value of the previous EM iterate and

Q(θ; θ (i) ) = -p(x 0:K |y 1:K , θ (i) ) log p(x 0:K , y 1:K |θ)dx 0:K . (2)
Then, the M-step consists in defining θ (i) as a minimizer for Q(•; θ (i) ), so that, by construction, the sequence (θ (i) ) i∈N decreases monotonically the ML cost function, and convergence guarantees to the ML estimate can be obtained under suitable assumptions. Note that, although the EM algorithm was originally an algorithm for computing ML estimates, it can also be modified for computation of MAP estimates, as we show in the next section.

PROPOSED GRAPHEM APPROACH

In this section we introduce the GraphEM algorithm for the estimation of the A matrix under sparsity constraints encoded in the prior distribution p(A) (see a discussion about the prior choice in Section 3.1). GraphEM aims at finding the maximum of p(A|y 1:K ) ∝ p(A)p(y 1:K |A), i.e., the MAP estimate of A. It is direct to show that the maximum of p(A|y 1:K ) coincides with the minimum of ϕ K (A) =log p(A)log p(y 1:K |A). In the model (1), ϕ K can be expressed recursively for k = 1, . . . , K:

ϕ k (A) = ϕ k-1 (A) -log p(y k |y 1:k-1 , A) (3)
= ϕ k-1 (A) + 1 2 log|2πS k (A)| + 1 2 z k (A) ⊤ S k (A) -1 z k (A), (4)
where z k (A) = y k -HAm k-1 (A) and S k (A) is the covariance matrix of the predictive distribution p(y k |y 1:k-1 ) = N (y k ; HAm k-1 (A), S k (A)), both being side products of the Kalman filter ran for a given A (see [START_REF] Sarkka | Bayesian Filtering and Smoothing[END_REF]Section 4.3]). Moreover, ϕ 0 (A) =log p(A) is the regularization function. The direct minimization of ϕ K is made difficult due to its recursive form. The EM approach allows to construct an upper bound for it, that is more tractable and thus easier to minimize. Let us denote

G k = P k (A ′ ) ⊤ (A ′ P k (A ′ ) ⊤ + Q) -1 the output of the RTS smoother for a given A ′ ∈ R Nx×Nx and set Σ = 1 K K k=1 P s k + m s k (m s k ) ⊤ , Φ = 1 K K k=1 P s k-1 +m s k-1 (m s k-1 ) ⊤ , C = 1 K K k=1 P s k G ⊤ k-1 + m s k (m s k-1 ) ⊤ .
Then, following [1, Section 12.2.3], we can show that the following function, parametrized by A ′ , majorizes the MAP objective function ϕ K for every A ∈ R Nx×Nx :

Q(A; A ′ ) = K 2 tr Q -1 (Σ -CA ⊤ -AC ⊤ + AΦA ⊤ ) + ϕ 0 (A) + C (5)
where tr is the trace operator, and C is a constant term independent from A. The GraphEM algorithm is then constructed as described in Alg. 1. The algorithm iterates alternating between an E-step and an M-step, ensuring the monotonic decrease of ϕ K (A (i) ) i∈N . In the E-step, the function

Q(A; A (i-1)
) is formed by running the Kalman filter and RTS smoother where the state matrix in the model is set to 1) . In the M-step, the matrix A (i) is updated by minimizing the latest generated majorizing function Q(A; A (i-1) ). In the following, we provide more details about the choice of the sparsity promoting prior and the M-step implementation.

A (i-

Algorithm 1 GraphEM algorithm

Initialization of A (0) For i = 1, 2, . . . (E-step) Run the Kalman filter and RTS smoother by setting A ′ := A (i-1) and construct Q(A; A (i-1) ) using Eq. ( 5).

(M-step) Update A (i) = argmin A Q(A; A (i-1) ) (see Section 3.2)

Choice of the prior

Sparsity is a key feature in statistical data processing to limit the degrees of freedom in parametric models. It has been widely used for graphical model inference [START_REF] Belilovsky | Learning to discover sparse graphical models[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical LASSO[END_REF], since a sparse matrix better allows to reveal interpretable and compact network of interdependencies between the entities. Since the ℓ 0 count measure is barely tractable, many statistical or learning approaches have invested in computable proxies, such as the ℓ 1 norm, also called Lasso penalty. In our problem, this amounts to set ϕ 0 (A) = γ A 1 with γ a positive penalty parameter. The ℓ 1 penalty is non-smooth, though convex and can be efficiently handled thanks to proximal optimization methods, as we will show hereafter.

Computation in the M-step

At an iteration i ∈ N, the M-step minimizes the function Q(A; A (i-1) ), obtained by plugging A ′ := A (i-1) in Eq.( 5). This minimization problem, sometimes called Lasso regression [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], has been much studied in the literature of optimization [START_REF] Schmidt | Fast optimization methods for l1 regularization: A comparative study and two new approaches[END_REF][START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF], and most of the methods proposed to solve it rely on the proximity operator [START_REF] Combettes | Proximal Splitting Methods in Signal Processing[END_REF][START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. 1 We recall that, for a proper, lower semi-continuous and convex function f : R Nx×Nx →] -∞, +∞], the proximity operator of f at A ∈ R Nx×Nx is defined as

prox f ( A) = argmin A f (A) + 1 2 A -A 2 F . (6)
Let us decompose Q(A;

A (i-1) ) = f 1 (A) + f 2 (A), where f 1 (A) = K 2 tr Q -1 (Σ -CA ⊤ -AC ⊤ + AΦA ⊤ )
, and f 2 = ϕ 0 is the prior that we have discussed in Section 3.1. Function f 1 is quadratic, so that, for θ > 0, we have

prox θf1 ( A) = mat Id ⊗ (KQ -1 ) + (θΦ -1 ) ⊗ Id -1 vec(KQ -1 CΦ -1 ) ,
due to the fact that the Lyapunov equation XA+AY = Z has for solution, A = mat((Id ⊗ X + Y ⊤ ⊗ Id) -1 vec(Z)), with vec and mat the re-ordering operations, using lexicographic order and ⊗ the Kronecker product. Note that if Q σ 2 Q Id, i.e., the state noise is isotropic, the expression simplifies into

prox θf1 ( A) = θK σ 2 Q C + A θK σ 2 Q Φ + Id -1
.

Second, the proximity operator for f is the simple soft thresholding operator:

prox θf2 ( A) = sign( A nm ) × max(0, | A nm |-θ) 1≤n,m≤Nx
Since function f 1 differentiable, a natural strategy would be to adopt the iterative thresholding method [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] or an accelerated version for it [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. Here, following the comparative analysis from [START_REF] Combettes | Proximal activation of smooth functions in splitting algorithms for convex image recovery[END_REF], we will prefer an algorithm that activates both terms via their proximity operator. More precisely, we will make use of the Douglas-Rachford (DR) algorithm, a fixed-point strategy for convex optimization that benefits from sounded convergence guarantees [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF] and has demonstrated its great practical performance in matrix optimization problems related to graphical inference applications [START_REF] Benfenati | A proximal approach for a class of matrix optimization problems[END_REF]. This leads to Alg. 2, that generates a sequence {A n } n∈N guaranteed to converge to a minimizer of f 1 + f 2 . In practice, for the iteration i in GraphEM, we run DR method with θ = 1 and initialization Z 0 = A (i-1) , i.e. the majorant function tangency point. Moreover, we stop the DR loop as soon as

|(f 1 + f 2 )(A n+1 ) -(f 1 + f 2 )(A n )|≤ ε (typically, ε = 10 -3 ). Algorithm 2 Douglas-Rachford algorithm for M-step Set Z 0 ∈ R Nx×Nx and θ ∈ (0, 2) For n = 1, 2, . . . A n = prox θf2 (Z n ) V n = prox θf1 (2A n -Z n ) Z n+1 = Z n + θ(V n -A n )

EXPERIMENTAL EVALUATION

The synthetic data are generated following model [START_REF] Sarkka | Bayesian Filtering and Smoothing[END_REF] with

N x = N y . A block-diagonal matrix A is considered, com- posed with b blocks of dimensions (b j ) 1≤j≤b , so that N x = b j=1 b j .
Each j-th diagonal block of A is randomly generated. Note that a projection on the space of matrices with maximum singular value less than one is applied, so that matrix A leads to a stable Markov process. Furthermore, we set H = Id, K = 10 3 and Q = σ 2 Q Id, R = σ 2 R Id, P 0 = σ 2 P Id with (σ Q , σ R , σ P ) some predefined values. Following this procedure, we created four datasets, whose parameters are provided in Table 1. Measures of quality to assess the performance of the methods are the relative mean square error (RMSE) on A, and the precision, recall, specificity, accuracy, and F1 score for the graph edge detection, when applying a threshold of 10 -10 on the entries of matrix A.

Dataset (b j ) 1≤j≤b (σ Q , σ R , σ P ) A ( 3

Table 1. Description of datasets

For each dataset, we ran GraphEM algorithm with the stopping rule |ϕ K (A (i) )ϕ K (A (i-1) )|≤ 10 -3 . In practice, a maximum number of 50 iterations is sufficient to reach this criterion. The algorithm is initialized setting a random and dense matrix A (0) , ensuring the stability of the process. Parameter γ, balancing the weight of the sparse prior, is optimized on a single realization thanks to a manual grid search, so as to maximize the accuracy score. We also provide the results obtained by the ML estimator, also computed using EM algorithm (MLEM). In this case, the M-step has a closed form solution [START_REF] Sarkka | Bayesian Filtering and Smoothing[END_REF]Th.12.5]. In addition, we compare with two Granger-causality approaches [START_REF] Bressler | Wiener-Granger causality: a well established methodology[END_REF] for graphical modeling. The first algorithm, pairwise Granger Causality (PGC) explores the N x (N x -1) possible dependencies among two nodes, at each time independently from the rest. The second approach, the conditional Granger Causality (CGC), operates similarly but for each pair of nodes, it also takes into account the information of the other N x -2 signals, in order to evaluate whether one node provides information to the other when the rest of signals are observed. We also do a manual grid search for finetuning the parameters of both PGC and CGC (more information can be found in [START_REF] Luengo | Hierarchical algorithms for causality retrieval in atrial fibrillation intracavitary electrograms[END_REF]). Note that PGC and CGC do not estimate a weighted graph but a binary one, so that RMSE is not calculated in those cases.

The results, averaged on 50 realizations, are presented in Table 2. MLEM does not promote sparsity in the graph which explains the poor results in terms of edge detectability. Moreover, GraphEM provides a better RMSE score on all examples. Regarding the graph structure, we can observe that GraphEM has also better detection scores, when compared with both PGC and CGC. We also display an example of graph reconstruction for dataset C in Fig. 1, illustrating the ability of GraphEM to recover the graph shape and weights.

CONCLUSION

This work proposes the GraphEM algorithm, an expectationminimization method for the estimation of the linear operator that encodes the hidden state relationships in a linear-Gaussian state-space model. Due to the use of a Lasso penalization term, we incorporate a sparsity constraint in our model. This is particularly well suited for modeling and representing the state entries interactions as a compact and interpretable graph. A proximal splitting algorithm is employed for solving the inner minimization problem of our EM approach. Numerical results illustrate the great performance of the method, when compared to several techniques for graphical modeling.

, 3, 3 ) ( 10 - 1 , 10 - 1 , 10

 310110110

Fig. 1 .

 1 Fig. 1. True graph (left) and GraphEM estimate (right) for dataset C.

Table 2 .

 2 Results for GraphEM, MLEM, PGC and CGC.

		method	RMSE	accur.	prec.	recall	spec.	F1
		GraphEM	0.081	0.9104	0.9880	0.7407	0.9952	0.8463
	A	MLEM PGC	0.149 -	0.3333 0.8765	0.3333 0.9474	1 0.6667	0 0.9815	0.5 0.7826
		CGC	-	0.8765	1	0.6293	1	0.7727
		GraphEM	0.082	0.9113	0.9914	0.7407	0.9967	0.8477
	B	MLEM PGC	0.148 -	0.3333 0.8889	0.3333 1	1 0.6667	0 1	0.5 0.8
		CGC	-	0.8889	1	0.6667	1	0.8
		GraphEM	0.120	0.9231	0.9401	0.77	0.9785	0.8427
	C	MLEM PGC	0.238 -	0.2656 0.9023	0.2656 0.9778	1 0.6471	0 0.9949	0.4198 0.7788
		CGC	-	0.8555	0.9697	0.4706	0.9949	0.6337
		GraphEM	0.121	0.9247	0.9601	0.7547	0.9862	0.8421
	D	MLEM PGC	0.239 -	0.2656 0.8906	0.2656 0.9	1 0.6618	0 0.9734	0.4198 0.7627
		CGC	-	0.8477	0.9394	0.4559	0.9894	0.6139

See also http://proximity-operator.net/


Maissa Sghaier

Emilie Chouzenoux

Jean-Christophe Pesquet

Serge Muller

A NEW SPATIALLY ADAPTIVE TV REGULARIZATION FOR DIGITAL BREAST TOMOSYNTHESIS

Keywords: Digital Breast Tomosynthesis, total variation, spatially adaptive regularization, 3D image reconstruction, optimization

Digital breast tomosynthesis images provide volumetric morphological information of the breast helping physicians to detect malign lesions. In this work, we propose a new spatially adaptive total variation (SATV) regularization function allowing to preserve adequately the shape of small objects such as microcalcifications while ensuring a high quality restoration of the background tissues. First, an original formulation for the weighted gradient field is introduced, that efficiently incorporates prior knowledge on the location of small objects. Then, we derive our SATV regularization, and integrate it in a novel 3D reconstruction approach for DBT. Experimental results carried out on both phantom and clinical data show the great interest of our method for the recovery of DBT volumes showing small lesions.

INTRODUCTION

In a will to reduce the mortality associated with breast cancer, heavy demands on the improvement of medical imaging modalities for breast cancer detection have been carried out in the last few years. Recently, Digital Breast Tomosynthesis (DBT) has demonstrated its superiority over X-ray mammography, the current modality of choice for breast cancer screening, for detecting lesions with comparable X-ray dose and a higher detection rate [START_REF] Pattacini | Digital mammography versus digital mammography plus tomosynthesis for breast cancer screening: the reggio emilia tomosynthesis randomized trial[END_REF].

One of the most important features of DBT is that it provides a quasi-3D image of the breast. This is of a crucial importance, especially for young women with heterogeneous breasts, from which cancer detection can be difficult when using standard mammography. In an analog way that it is a difficult task to see a needle in a haystack, detecting by visual inspection a cancer lesion in a standard 2D mammography remains quite challenging. Yet, three-dimensional DBT provides a piercing view through the haystack, i.e., the breast, by reducing the impact of tissue structures overlapping [START_REF] Sechopoulos | A review of breast tomosynthesis. Part I. The image acquisition process[END_REF] and so potentially increases the sensitivity and specificity of cancer detection. A key ingredient to guarantee such benefits is the image reconstruction algorithm. Exhaustive investigations and comparative studies on DBT reconstruction algorithms have clearly shown the advantages of iterative reconstruction algorithms [START_REF] Sechopoulos | A review of breast tomosynthesis. Part II. Image reconstruction, processing and analysis, and advanced applications[END_REF]. Many of these algorithms rely on a convex optimization formulation employing a total variation (TV) regularization function [START_REF] Sidky | Image reconstruction in digital breast tomosynthesis by total variation minimization[END_REF][START_REF] Sidky | Enhanced imaging of microcalcifications in digital breast tomosynthesis through improved imagereconstruction algorithms[END_REF][START_REF] Metin | An iterative tomosynthesis reconstruction using total variation combined with non-local means filtering[END_REF]. The TV regularizer has illustrated its great performance in various applications in image processing [START_REF] Condat | Discrete total variation: New definition and minimization[END_REF][START_REF] Abergel | The shannon total variation[END_REF][START_REF] Chambolle | An upwind finitedifference method for total variation-based image smoothing[END_REF]. However, the undesired staircasing side-effect of TV [START_REF] Nikolova | Local strong homogeneity of a regularized estimator[END_REF] has led researchers to propose alternative TV-based penalties, by proposing various smooth approximations [START_REF] Aujol | Some first-order algorithms for total variation based image restoration[END_REF][START_REF] Blomgren | Total variation image restoration: numerical methods and extensions[END_REF][START_REF] Chan | High-order total variation-based image restoration[END_REF][START_REF] Chan | Image decomposition combining staircase reduction and texture extraction[END_REF], by investigating the non-locality [START_REF] Chierchia | A nonlocal structure tensor-based approach for multicomponent image recovery problems[END_REF], or by computing a posterior mean instead of the MAP estimation framework [START_REF] Louchet | Total variation denoising using posterior expectation[END_REF]. Better recovery of anisotropic objects can also be reached, using anisotropic TV-based regularization [START_REF] Luo | Adaptive weighted total variation minimization based alternating direction method of multipliers for limited angle ct reconstruction[END_REF][START_REF] Merveille | A variational model for thin structure segmentation based on a directional regularization[END_REF][START_REF] Liao | An image denoising fast algorithm for weighted total variation[END_REF][START_REF] Xu | Weighted total variation based convex clustering[END_REF][START_REF] Biton | Adaptive anisotropic total variation-a nonlinear spectral analysis[END_REF]. Nonetheless, to the best of our knowledge none of them have provided a TV formulation that takes into consideration the heterogeneity (in shape and texture) within the image to be regularized. Breasts consist of different anatomical components such as tissues (background) and lesions that require a more sophisticated regularization when reconstructing the DBT volume. In such challenging context, our contribution is to propose a new Spatially Adaptive TV (SATV) regularization function which responds to the need for heterogeneous spatial regularization in DBT reconstruction, and accounts for the different morphological contents of the breast. Henceforth, we provide a new definition of the gradient field in the discrete image that takes into consideration an input detection map incorporating prior knowledge on the structural contents of the image. Thus, the SATV operator acts differently according to the sought local information around each voxel.

The paper is organized as follows: Section 2 provides the description of our convex optimization framework for DBT reconstruction. In Section 3, we present our main contribution which lies in the mathematical formulation of the SATV regularization function, and we discuss its incorporation in 3D reconstruction algorithms. Section 4 shows the experimental results carried out on both physical phantom and clinical data and discusses the obtained qualitative results. Finally, conclusions are drawn in Section 5.

CLINICAL TASK-BASED RECONSTRUCTION IN DBT

From a mathematical point of view, the reconstruction of the volumic image from DBT measurements is an inverse problem. The forward model is often simplified into the following linear form [START_REF] Sechopoulos | A review of breast tomosynthesis. Part II. Image reconstruction, processing and analysis, and advanced applications[END_REF]:

p = Ad + e (1)
where p ∈ R n is a vector derived from the acquired projections, d ∈ R m is a vector representing the unknown volume, A ∈ R n×m is a matrix describing the geometry of the system, and e ∈ R n is the noise. The objective is then to have an accurate reconstruction d from under-sampled data p. Since the sampling rate for acquiring the projections p is very low in standard DBT configurations, it is challenging to obtain a solution d which reaches a globally satisfactory visual quality. It appears also crucial to take into consideration, in the reconstruction process, the different clinical tasks that are carried out by the radiologist when reading DBT images, in particular microcalcifications detection. To reach both of these goals, we formulate the solution to (1) as the minimization of a penalized least squares cost function:

minimize d∈[0,dmax] m 1 2 p -Ad 2 2 + Φ(d) -D(d). (2)
Hereabove, Φ : R m → R is a TV-based spatial regularizer, and D : R m → R denotes a detectability function enhancing the visibility of certain objects related to the sought clinical task. In a previously published work we have shown that the minimization of the cost function (2) yields good results in terms of microcalcifications enhancement in DBT reconstruction [START_REF] Sghaier | A new approach for microcalcification enhancement in digital breast tomosynthesis reconstruction[END_REF] while a slight spreading on the microcalcifications shapes can be observed. Furthermore, it can be expected that a better image quality trade-off is obtained by designing a TVbased regularization taking into account some prior knowledge on the different regions of the volume, as it is already done for building the clinical-based term D. Subsequently, we introduce a new total variation operator that enables an apter spatial regularization in the DBT volume accounting for the sought morphological content in the neighborhood of each voxel.

SATV REGULARIZATION

Mathematical formulation

With the aim to make regularization adaptive to the local contents of the DBT volume d ∈ R m , we define the new SATV regularization function as follows

Φ(d) = m i=1 ψ (1 -η)λ i (∆d) i + η(1 -λ i )δd i , (3)
where δ = [1 1 1] , ψ : R 3 → R represents a sparsity promoting function which applies either on the image intensity or its gradients, η controls the overall contribution of each regularizing term, ∆ denotes the 3D gradient operator that is expressed as

(∆ i ) = [∆ x i ∆ y i ∆ z i ] (4)
where ∆ x i (resp. (∆ y i ) and (∆ z i )) denotes the horizontal, (resp. vertical and depth) gradient at voxel i ∈ {1, . . . , m}, with zero-boundaries assumption. It must be emphasized that the nonnegative weight parameters (λ i ) 1≤i≤m play a crucial role in the adaptation to the local spatial contents which is performed. In particular when, for every i ∈ {1, . . . , m}, λ i = 0 and ψ = • 1 , Φ reduces to the 1 regularization, whereas, when for every i ∈ {1, . . . , m}, λ i = 1 and ψ = • , we recover the standard isotropic total variation. In the case of DBT volumes, we would like to apply different regularizations in the areas where microcalcification are present/absent. The advantage of formulation ( 2) is that it allows us to exploit the detection map, extracted from a Computer-Aided Detection like tool [START_REF] Christoyianni | Computer aided diagnosis of breast cancer in digitized mammograms[END_REF], that enables the retrieval of the location of microcalcifications in the DBT volume. Using this detection map, a prior estimate µ i of the normalized intensity value of voxel i ∈ {1, . . . , m} can be derived [START_REF] Sghaier | A new approach for microcalcification enhancement in digital breast tomosynthesis reconstruction[END_REF]. We recall that the microcalcifications are defined, in this framework, as localized structures lying within a predefined radius range. Henceforth, the proposed expression of the weight parameters is given by

λ i =      1 if µ i ≤ ν -θ θ + ν -µ i 2θ if ν -θ ≤ µ i ≤ ν + θ 0 otherwise, where ν ∈ [0, +∞[ with θ ∈]0, min(ν, 1 -ν)].
The value µ i is high in zones where structures are present, while it is low in other areas. The threshold ν allows us to set the binary transition between the two zones, while the parameter θ is used to smooth the latter transition.

Reconstruction algorithms

In this section, we briefly explain how Problem (2) can be numerically handled when SATV regularization is used. From (3), function Φ can be reexpressed as

∀d ∈ R m , Φ(d) = Ψ (Ld), (5)
where

∀u = (u i ) 1≤i≤m ∈ (R 3 ) m , Ψ (u) = m i=1 ψ u i ), (6)
and

L = (1 -η)(Λ ⊗ δ)∆ + η(1 -Λ) ⊗ δ ∈ R 3m×m , ( 7
)
where ⊗ denotes the Kronecker product and Λ = Diag(λ 1 , . . . , λ m ). The function ψ can be assumed proper, lower-semicontinuous and convex and the detectability function D can be chosen linear [START_REF] Sghaier | A new approach for microcalcification enhancement in digital breast tomosynthesis reconstruction[END_REF]. Since Φ is the composition of matrix L and the proper, lower-semicontinuous and convex function Ψ , a variety of splitting iterative algorithms can be applied to solve Problem (2):

• If the proximity operator of ψ has a closed form expression, proximal algorithms [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] which include ADMM [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] and primal-dual methods [START_REF] Komodakis | Playing with duality: An overview of recent primal? dual approaches for solving largescale optimization problems[END_REF] provide efficient solutions.

• If ψ is Lipschitz-differentiable, then the problem can even be solved in a simpler manner through (possibly accelerated) projected gradient algorithms, where a projection onto the hypercube [0, d max ] m is performed at each iteration.

EXPERIMENTAL RESULTS

In this section, we present experimental results carried out first on a physical phantom, then on two realistic clinical datasets. Each dataset contains 9 projections acquired using a DBT commercial system (Senographe Essential, GE Healthcare) with an angular aperture of 25 o . The detector is composed of 3062 × 2394 detector elements of 100µm 2 in size. The volumes are reconstructed under a 100µm × 100µm × 1mm sampling grid. We aim at identifying small microcalcifications of size within the range [0.2, 0.4]mm from the rest of the breast. Therefore, we rely on a map µ allowing us to detect structures with the latter size range following the approach proposed in [START_REF] Sghaier | A new approach for microcalcification enhancement in digital breast tomosynthesis reconstruction[END_REF]. We then construct (λ i ) 1≤i≤m as detailed in Section 3.1. The detection map is binarized by setting the threshold value ν so that the geometric shape of the detected microcalcifications are well recovered. Moreover, θ is manually finetuned to optimize the transitions between breast background and microcalcifications. Finally, we set η in order to achieve a certain image quality optimization of DBT. In the general expression (3) of SATV regularization, a smooth approximation of the Euclidean norm is used as a penalty function, namely for every

u ∈ R 3 , ψ(u) = β u 2 + 2
, where (β, ) ∈ (R + ) 2 control the weight and smoothness of the regularization, respectively. The upper bound in the range value constraint is chosen equal to d max = 4095. As explained in Section 3.2 the minimization of (2), in the considered setting, can be performed by projected gradient descent. The algorithm was initialized with the Filtered Back Projection solution, and ran until a maximum number of iterations (here equal to 200), which was shown to be sufficient to reach practical convergence in our experiments.

Physical phantom data

We have first evaluated our approach using the physical phantom accreditation Model 015 (CIRS). Its dimension is 45 × 102 × 108mm 3 , with uniform background. For the sake of visualization, we only display a zoomed region containing 6 microcalcifications of 0.3mm in size and with known locations. We display in Fig. 1(a) the constructed (λ i ) 1≤i≤m as explained in Section 3.1. We set ν = 1.5 × 10 -3 and θ = ν/3. It is worth mentioning that the background is merely equal to zero in the detection map. Yet, the detected microcalcifications present different grey level values which will lead to a different regularization magnitude for each of them. The benefits of our proposed method as illustrated in Fig. 1(d) is twofold. First, compared to the results obtained with a classical DBT reconstruction (Fig. 1(b)), we significantly enhance the visibility as well as the sharpness of the detected microcalcifications. Second, compared to the results obtained with a clinical task-based reconstruction using a classical TV, as illustrated in Fig. 1(c), we greatly improve the shape restoration of the enhanced microcalcifications. We emphasize that the comparison between the SATV-based reconstruction and the classical TV-based reconstruction is fair, since the regularization parameters, namely β and η, are set such that the background is restored with the same quality.

In this example, β = 600, = 10 for classical TV while β = 630, = 10, η = 0.048 in the case of the new SATV regularization.

Clinical data

We have also assessed the performance of our proposed approach on a clinical scenario. We considered here two BI-RADS b breast density clinical datasets, which correspond to breasts with scattered areas of fibroglandular density [START_REF]Acr bi-rads atlas: breast imaging reporting and data system[END_REF]. We show in Fig. 2 a case compromising a potentially malignant cluster of microcalcifications, while in Fig. 3 a case containing an isolated microcalcification. We display the map of (λ i ) 1≤i≤m obtained from the estimated signal where ν = 0.09 and θ = ν/3, in Fig. 2(a) and Fig. 3(a) respectively. Despite the presence of false positives in these latter maps, we note that the enhancement is more effective in the lesion area. One can clearly observe in Fig. 2 and in Fig. 3 the behavior already unveiled in the physical phantom dataset. With our proposed method, we highlight the advantage of simultaneously enabling a robust preservation of the shape of the enhanced microcalcifications and of a high quality restoration of the background tissues. Hereagain, the reconstruction with our proposed approach and the one based on classical TV provide the same breast texture quality, after fine-tuning of the regularization parameters η and β for this purpose. Here, β = 600, = 10 for classical TV and β = 660, = 10, η = 0.09 when SATV penalization is used.

CONCLUSION

Despite the tremendous amount of research works highlighting the success of TV in various image processing applications, this penalty does not model properly the heterogeneity within images. In this paper, our main contribution lies in a novel definition of a spatially adaptive total variation which incorporates prior shape knowledge on the image to be restored. In the DBT context, the proposed SATV allows to take into account the expected location of small microcalcifications. We have demonstrated on different datasets the interest of our proposed approach with respect to classical TVbased DBT reconstruction. One future direction for further improvements is to investigate techniques for a more automatic setting of the different parameters involved in SATV, depending on the targeted task and images to be processed.

Fig. 1 :

 1 Fig. 1: Region of ACR phantom slice containing 6 ROIs: (a) Weights (λi) 1≤i≤m , (b) DBT reconstruction with classical TV and no clinical term, (c) DBT reconstruction with classical TV and clinical term (d) DBT reconstruction with SATV regularization and clinical term.

Fig. 2 :

 2 Fig. 2: Sample of slices of a DBT reconstruction with size 3062 × 994 × 67(pixels) : (a) Constructed (λi) 1≤i≤m from the representative clinical data. (b) DBT reconstruction with classical TV regularization and without the clinical term, (c) DBT reconstruction with classical TV regularization and the clinical term, (d) DBT reconstruction with SATV regularization and the clinical term.

Fig. 3 :

 3 Fig. 3: Sample of slices of a DBT reconstruction with size 2344 × 868 × 44(pixels) : (a) Constructed (λi) 1≤i≤m from the representative clinical data. (b) DBT reconstruction with classical TV regularization and without the clinical term, (c) DBT reconstruction with classical TV regularization and the clinical term, (d) DBT reconstruction with SATV regularization and the clinical term.


Danilo Longo
email: danilo.longo@synchrotron-soleil.fr

Marie-Laure Bocquet

Nicolás Lorente

Hervé Cruguel

François Debontridder

Sébastien Royer

Pascal David

Alexandra Palacio-Morales

Tristan Cren

Nadine Witkowski

Christophe Brun

¶pasteur

Surface Phase Nucleation of Lead Monoatomic Layers on Si(111) Induced by Manganese Phthalocyanine Molecules

Keywords: MPc=metal-phthalocyanine surface phase transition, metal/semiconductor interface, lead monolayer, molecule-substrate interaction, phthalocyanines, two-dimensional metallic systems, molecular electronic/magnetic properties

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The recent search for topological superconductivity 1 prompted intense eorts to study the basic coupling between magnetic impurities and superconducting materials. 25 Recently, it has been shown that this coupling strongly depends on the dimensionality of the superconductor 6 and that 2D materials might be much more promising than 3D ones for such purpose.

79 In this regard, metallic monolayers of either Pb or In epitaxially grown on semi- conducting substrates constitute a very interesting class of systems.

1012 Moreover, coupling these 2D materials with arrays of metallo-organic molecules, e.g. metal-phthalocyanines (MPcs), seems an appealing way to realize hybrid magnetic/superconducting systems.

One of the advantages of using MPcs is the possibility of tailoring the molecule-substrate coupling by selecting the most appropriate metallic atom M to coordinate within the molecular structure. 13 For this reason the interaction of MPcs with dierent substrates has been extensively investigated.

1419 Interestingly, the set of these studies reveals that, in addition to the coordinated atom M, the substrate also plays a crucial role in balancing the relative emphasize the rich and new behavior of hybrid organic/metal systems emerging when the dimensionality of the metallic substrate is reduced. In particular, this could enable a new way of tuning the coupling between magnetic molecules and a 2D superconductor, with the aim of stabilizing well-dened magnetic ordering of the molecular array.

Results and discussion

Experimental results

The adsorption of MnPc molecules on Pb monolayers was carried out after proper preparation of the substrate. Single atomic layers of lead were obtained by rst depositing 1.65 monolayer (ML) of Pb on top of a freshly annealed Si(111) surface with 7 × 7 reconstruction. Then, the sample was annealed at about 600 K for a few minutes. Finally, the crystalline quality of the Pb/Si(111) surface was checked with both STM and LEED (see Methods section for details).

Figure 1a shows STM and LEED results of the 3033 Accordingly, the complex LEED pattern of Figure 1a is due to the overlap of the diraction spots arising from the Si(111)-1 × 1 substrate and the three √ 7 × √ 3-Pb domains (one of them is highlighted by the purple arrows). [START_REF] Stepanovsky | The Dense α-√ 3 × √ 3Pb/Si(111) Phase: A Comprehensive STM and SPA-LEED Study of Ordering, Phase Transitions and Interactions[END_REF][START_REF] Ganz | Growth and Morphology of Pb on Si( 111)[END_REF][START_REF] Choi | High-Resolution Core-Level Photoemission Study of Dense Pb Overlayers on Si(111)[END_REF] The STM image of gure 1a shows only one of these

Si-1×1

√7×√3-Pb

Si-1×1

√3×√3-Pb

Si-1×1

√3×√3-Pb

Si-1×1

b) SIC c) MnPcs-induced SIC a) √7×√3-Pb [START_REF] Horikoshi | Structural Phase Transitions of Pbadsorbed Si(111) Surfaces at Low Temperatures[END_REF][START_REF] Stepanovsky | The Dense α-√ 3 × √ 3Pb/Si(111) Phase: A Comprehensive STM and SPA-LEED Study of Ordering, Phase Transitions and Interactions[END_REF]3740 According to this, the LEED pattern of Figure 1b shows the image (V T =+0.5 V and I T =50 pA) in (e) shows MnPcs self-assembled in a trimer. The relative orientation of a timer with respect to the substrate is depicted in the drawing (f ).

√ 3 × √ 3-
All STM images were acquired at 77 K at which temperature tip induced molecular motion exists.

spectively. A clear MnPcs-induced SIC reconstruction is observed in both cases. Moreover, one sees that isolated molecules and self-assembled trimers are the most common congurations observed at low molecular densities. Remarkably, this self-assembly behavior of MnPc molecules on the Pb monolayer seem to be independent of the substrate temperature in the range 100 < T < 300 K. This observation suggests that molecules do diuse during adsorption and that their nal conformations are not limited by diusion barriers. Our study also reveals that, at low molecular coverage, MnPcs do not self-assemble to form extended long-range ordered molecular arrays on Pb/Si(111) monolayers. This is in contrast to the spontaneous assembly of squared lattices of MnPcs usually observed on Pb thin lms or bulk monocrystals.

1419,43 Moreover, to our knowledge, trimer arrangement of MnPcs has not been previously reported.

1416,18,19 Furthermore, in contrast to what has been previously reported in the case of C 60 molecules deposited on the Pb/Si(111) monolayer, [START_REF] Matetskiy | Structural Transformations in Pb/Si(111) Phases Induced by C 60 Adsorption[END_REF] MnPc molecules lie on defect-free regions.

Not only do the molecules remain isolated or arranged in trimers, but they present precise orientations with respect to the substrate. Similar molecular organization was also found on the pristine SIC surface. Figure 3 shows ∼1.3 molecules/100 nm 2 of MnPcs deposited on the SIC phase kept at room temperature.

Here the SIC surface structure is not aected by MnPcs adsorption. However, MnPcs retain the same self-assembly properties than the ones observed on the MnPcs-induced SIC phase, that is single MnPcs and self-assembled trimers are the only molecular congurations observed on the Pb monolayer surface (see inset of gure 3) and their arrangement reects the threefold symmetry of the Pb substrate.

Computational results

To elucidate the details of the macroscopic 5053 However, after adsorption on the Adsorption of MnPcs on Pb: 3D vs 2D substrate

√ 7× √ 3-
In order to get more insight on the link between substrate dimensionality and moleculesubstrate coupling strength, we compare below our results with those already reported for adsorption of MPc molecules on 3D substrates. Adsorption of MnPcs on Pb(111) has already been extensively investigated. However, none of these studies has ever reported anything similar to what we observed on the Pb monolayer.

MnPcs deposited on Pb(111) self-assemble in close-packed square lattices even at submonolayer coverage suggesting that their self-assembly is essentially driven by the moleculemolecule interaction as well as the fourfold symmetry of the molecules themselves.

1416,43

The weak molecule-substrate coupling manifests itself by the emergence of dierent superstructures superimposed on the square lattice of MnPcs. At present dierent computational approaches exist to describe the adsorption of a MnPc molecule on bulk Pb(111).

5053 All studies agree that the MnPc molecule in the gas phase has an electronic conguration that gives rise to a total spin S=3/2. After adsorption on Pb(111), the magnetic orbital the most strongly hybridized with the substrate is the Mn d z 2

and any charge transfer always occurs from the Pb(111) substrate to the MnPc molecule.

This is in agreement with our ndings. Then, depending on the microscopic details included in the model, e.g. electron correlations or hybridization, dierent predictions on the charge transfer and magnetic state of MnPcs after adsorption come out. In agreement with their new STS results, Minamitani et al.

53 reported about a collective magnetic state delocalized on the organic ligand and not localized only on the central Mn atom. According to this, they have shown that after adsorption on Pb(111) the total molecular spin decreases from S=3/2 to S=1 due to a charge transfer of 2.33 electrons essentially localized on the organic ligand and a small contribution of the Mn atom with 0.1 electrons. Our DFT analysis agrees with this result and reveals that the magnetism arising from the d z 2 orbital survives after adsorption of MnPcs on the √ 7 × √ 3-Pb monolayer despite its strong hybridization with the substrate while the magnetism arising from d xz /d yz orbitals (that are coupled to the organic ligands) is suppressed by the charge transfer. This result is dierent from those reported by previous studies in which the total spin reduction to S=1 or even S=1/2 was interpreted as an interaction with substrate/ligands strictly localized on the magnetic Mn atom and not involving organic ligands. However, the molecular-induced surface reconstruction of Au(111) occurs at a molecular coverage of 1 ML. This is in striking contrast with our experiment in which we observed a long-range surface reconstruction of Pb monolayer driven by a dilute amount of MnPcs. We believe that this long-range eect relies on the fact that the Pb monolayer is a particular type of strained surface. Indeed, it was shown that long-range repulsive interactions within the Pb monolayer are responsible for the emergence of a huge number of structural phases.

59

These repulsive interactions are also responsible for macroscopic rearrangement of phases with dierent symmetries when the Pb coverage is slightly changed. This fact suggests that a local perturbation (even small) of the Pb monolayer may actually be felt on macroscopic distances through variations of the surface strain inside the system.

Another important eect that is well known to occur in proximity of an adsorbed molecule is the charge transfer which induces a redistribution of the surface electrons within the substrate.

58,60,61 A direct consequence of this electronic reorganization is a variation of both surface strain and local atomic arrangement at the molecule-substrate interface. Interest-ingly, our DFT calculations show that a strong charge transfer takes place from the Pb monolayer to the MnPc molecule. Furthermore, our charge density maps reveal that the transferred charge is mostly accumulated at the Mn-Pb interface (see gure S2 in supporting information).

Taking into account the various features of our surface system we tentatively propose a mechanism behind the macroscopic phase transition of the Pb monolayer induced by MnPcs adsorption. We believe that the charge transfer between Pb and MnPcs may act as local perturbations that aect the surface strain of the whole system and initiate the structural transformation. The We have also shown that there exists intrinsic dierences between the electronic/structural properties of the Pb monolayer and those of thicker Pb substrates that might be ascribed to the participation of the Pb/Si(111) interface in the molecule-substrate interaction. In particular, a dierent self-assemblying behavior as well as a dierent electronic/magnetic state of individual MnPcs was found in the present study. This interesting result opens up the possibility of using the substrate dimensionality to tune the molecule-substrate coupling as well as the electronic and magnetic properties of organic hybrid systems.

Our results can have a strong impact in the search for organized magnetic nanostructures on superconductors. In particular, the rich and new behavior of hybrid organic/metal systems occurring when the dimensionality of the metal substrate is reduced enables a new way of tuning the coupling between single and assembled magnetic molecules and a 2D superconductor. This could be important for research aiming to stabilizing magnetic order in metal-organic systems.

Methods

Sample preparation

All the experiments were carried out in Ultra-High Vacuum (UHV) with a base pressure bet-

7 × √ 3 - 7 × √ 3 - 7 × √ 3 - 7 × √ 3 - 7 × √ 3 -

 7373737373 strength of the molecule-molecule and molecule-substrate interaction. As an example, it has recently been shown that MPcs deposited on the √ In phase of the In/Si(111) bilayer form closely packed square lattices 19 indicating the prevalence of the molecule-molecule inter- action. As small as it is, the molecule-substrate interaction remains non-negligible since the critical superconducting temperature (T C ) of the In bilayer increases or decreases according to the particular MPc (M=Mn or Cu) deposited on the surface. If self-assembly of MPcs on the In bilayer results in square lattices, the adsorption of molecules that interact more strongly with the substrate might lead to very dierent behavior of the hybrid system. For instance, the deposition of 7,7,8,8-tetracyanoquinodimethane (TCNQ) molecules on top of the √ 7 × √ 3-In/Si(111) bilayer leads to the formation of a surface metal-organic network (SMON) with inclusion of In atoms extracted from the underlying metallic substrate. 20 Simultaneously, the metallic area surrounding the SMON undergoes a surface reconstruction from the √ In phase to a novel 5 × 5-In phase as a consequence of the local decrease in In coverage. Surface reconstruction of a 2D metallic system induced by molecular adsorption was also reported by Matetskiy et al. 21 In that study, C 60 molecules were deposited on top of the √ Pb phase of the Pb/Si(111) monolayer. Self-assembly of disordered C 60 islands was accompanied by the expulsion of Pb atoms from the monolayer allowing the molecular domains to be incorporated within the layer. Expelled Pb atoms were subsequently reintegrated in the regions surrounding the C 60 islands by increasing the local Pb coverage thus leading to a structural phase transition of the Pb surface from the √ Pb to the denser hexagonal incommensurate HIC phase. Here, we report on a combined Scanning Tunneling Microscopy (STM), Low Energy Electron Diraction (LEED) and Density Functional Theory (DFT) study on the properties of crystalline Pb monolayers grown on Si(111) in presence of minute concentrations of manganese phthalocyanines (MnPcs). We present a new structural phase transition of the Pb/Si(111) monolayer, from the √ 7× √ 3-Pb/Si(111) to the Striped Incommensurate (SIC) phase, taking place upon adsorption of ∼0.18 molecules/100 nm 2 of MnPcs. In contrast to previous studies on similar systems, we found that the surface reconstruction mechanism does not imply the expulsion/reintegration of Pb atoms close to the molecules but is rather related to the strong local molecule-substrate interaction. The structural transition is ultimately promoted by the strained nature of the Pb phases and the energetic stability of the MnPc/SIC-Pb/Si(111) system compared to the MnPc/ √ Pb/Si(111). Moreover, LEED measurements reveal that ∼0.18 molecules/100 nm 2 of MnPcs are already able to nucleate the new phase. This molecule-induced surface reconstruction as well as the observation of a frustrated molecular self-assembly reecting the threefold symmetry of the substrate reveal an enhanced moleculesubstrate coupling stronger than what is known on bulk or ultrathin Pb substrates. This is corroborated by our DFT calculations showing that the Pb/Si(111) interface plays a crucial role in determining the molecule-substrate coupling strength as well as the physical properties of the system such as the magnetic state of the MnPc molecules. Ultimately, our results

3 - 3 -

 33 Pb domains and reveals the typical periodic arrangement of Pb atoms in straight lines. Pb unit cell contains 6 Pb atoms per 5 Si atoms corresponding to a Pb coverage of 6 5 ML=1.2 ML.

Figure 1 : 3 -× √ 3 -

 133 Figure 1: (a) and (b) show STM images (top) and LEED patterns (bottom) of a Pb/Si(111) monolayer prepared in the √ 7× √ 3-Pb and SIC phase respectively. (c) shows the STM image (top) and the LEED pattern (bottom) of the SIC reconstruction induced by adsorption of ∼0.18 molecules/100 nm 2 of MnPcs on the Pb monolayer initially prepared in the √ 7 × √ 3-Pb phase (substrate kept at room temperature during MnPcs deposition). All STM images were measured at 77 K over a 50×50 nm 2 area: (a) V T =+1 V / I T =20 pA, (b) V T =-50 mV / I T =80 pA and (c) V T =+0.1 V / I T =40 pA. Corresponding LEED patterns recorded at E∼80 eV and (a) 100 K (b) 300 K (c) 100 K.

3 - 7 × √ 3 -Figure 2 : 7 × √ 3 -

 373273 Figure 2: 46×46 nm 2 STM images in (a) and (b) (V T =+0.3 V and I T =20 pA and V T =+0.4 V and I T =80 pA respectively) show ∼0.7 molecules/100 nm 2 and ∼4.6 molecules/100 nm 2 of MnPcs respectively adsorbed on the MnPcs-induced SIC surface. MnPcs were deposited on the √ 7 × √ 3-Pb surface kept at approximately 150 K. The 5×5 nm 2 high-resolution STM image (V T =+0.4 V and I T =80 pA) in (c) shows a single MnPc molecule corresponding to the conguration I in the sketch (d). The schematic drawing in (d), deduced from the experimental data, represents the three adsorption orientations of single MnPcs with respect to the Pb surface (purple arrows). Red arrows indicate the two nitrogen atoms within MnPcs through which a substrate direction passes. Black segments highlight one of the main molecular axis of a MnPc in a given conguration. The 6×6 nm 2 high-resolution STM

  Figure 2c displays an isolated MnPc molecule adsorbed on the Pb monolayer surface. The sketch of gure 2d, resulting from statistical analysis of several STM images, shows that isolated MnPcs can be found in only three possible orientations, rotated by an angle of 120 • between each other, with respect to the symmetry directions of the Pb monolayer (indicated by purple arrows). The MnPc molecule in gure 2c, for instance, corresponds to the conguration I depicted in gure 2d. For each of these molecular congurations one of the substrate directions is always aligned with the axis passing through the two nitrogen atoms indicated by the red arrows (see gure 2d) within the experimental error of about 2 • . All these results and considerations indicate a strong and local inuence of the Pb monolayer on the adsorption conguration of MnPcs. This is in stark contrast with the lack of surface registry found for MnPcs adsorption on the (111) surface of bulk Pb 15,16,43 or Pb thin lms. 14 The sketch of gure 2f, deduced from a statistical analysis of STM images, depicts the arrangement of MnPcs inside self-assembled trimers. The registry of the trimer with respect to the Pb surface high-symmetry directions is indicated by purple arrows. One measured trimer is presented in gure 2e. Comparison between gure 2d and gure 2f reveals that the single MnPcs inside a trimer always retain the same orientation with respect to the Pb substrate as do isolated molecules. More interestingly, all trimers observed in our experiments always consisted of MnPc molecules with the three possible orientations depicted in gure 2d.

Figure

  Figure 2f also shows that each Pb direction passes through Mn atoms of neighboring MnPcs within the experimental error of about 2 • . Moreover, these Mn atoms constitute the vertices

  Figure 3: 60×60 nm 2 STM image (V T =-50 mV and I T =50 pA) showing ∼1.3 molecules/100 nm 2 of MnPcs deposited on top of the pristine SIC phase of the Pb/Si(111) monolayer kept at room temperature. MnPcs remain isolated or self-assemble in trimers (see inset: 7×7 nm 2 , V T =-50 mV and I T =50 pA). Both images were measured at 300 mK.

Figure 3 3 -

 33 also show a dimer of MnPcs. These congurations were rarely observed during the experiments and they can be seen as trimers in which a molecule is missing. The similar behavior of MnPcs on the SIC and the MnPcs-induced SIC phase suggests that the molecule-substrate interactions are similar on both surfaces. This observation further conrms that the MnPcs adsorption on top of the Pb monolayer actually transforms the √ 7 × √ 3-Pb phase into the SIC phase. In summary, at low coverage MnPc molecules deposited on the √ 7× √ Pb monolayer induce a surface reconstruction into the incommensurate SIC phase. The above data show that molecules form strong and local bonds with the substrate such that even if they diuse during deposition, the molecules do not nd the required adsorption sites on the incommensurate phase to form large and ordered assemblies. This forces the molecules to lie isolated or in trimers and always in congurations that reect the threefold symmetry of the Pb surface.

3 × √ 3 -

 33 structural phase transition of the Pb monolayer induced by MnPc molecules adsorption we performed DFT calculations. Here, one should keep in mind that the SIC phase is an incommensurate surface reconstruction and it cannot be exactly treated with electronic-structure calculations. This is the reason why we assumed a commensurate √ Pb atomic arrangement instead of the actual incommensurate SIC phase. Note that this simplication is commonly done to address the SIC electronic properties with DFT methods. 10,26,39,44 We did two separate calculations, rst the adsorption of a single MnPc molecule onto the √ 7 × √ 3 phase and second, onto the √ 3 × √ 3 phase. The minimum-energy conformations of the two systems are displayed in gure 4.

Figure 4 : 3 - 7 × √ 3 - 7 × √ 3 - 3 × √ 3 - 3 - 3 × √ 3 - 2 Figure 5 : 3 -

 43737333333253 Figure 4: DFT optimized structure and Terso-Hamann STM topography 4547 at xed

3 → 7 × √ 3 - 3 × √ 3 - 7 × √ 3 - 3 -

 37333733 SIC structural phase transition of the Pb monolayer. MnPc molecules do stabilize better on a denser Pb phase and favors the reconstruction into an incommensurate phase. The initial Pb arrangement of the √ Pb phase proceeds towards a reorganization to locally form denser areas of √ Pb and other areas with Pb troughs leading to the overall incommensurate SIC surface structure. In addition, the well-known fact that Pb adatoms in the overlayer are soft and easily mobile 37,39 might help facilitating this process.DiscussionThe present LEED/STM study reveals a new long-range eect taking place when MnPc molecules are adsorbed at low coverage on a Pb/Si(111) monolayer. In particular, we have shown that the adsorption of a very dilute amount of MnPcs is sucient for Pb atoms initially arranged in the √ Pb phase to reorganize themselves into the SIC phase on a macroscopic scale. On this MnPcs-induced SIC surface, the molecules remain either isolated or self-assembled in trimers and their adsorption conguration reects the threefold symmetry of the substrate. Moreover, our DFT calculations have shown that both the energy gain and the charge transfer of single MnPc molecules adsorbed on a pure Pb phase supporting the experimentally observed structural transition.

7 × √ 3 -

 73 53 of a single MnPc molecule adsorbed on 3D Pb(111), show that the spin density remains distributed around the molecule and does not penetrate the Pb substrate, in contrast to the present calculations.We would also like to discuss the interesting case of Au(111), Ag(111) and Cu(111) 3D substrates, all having 2D surface states at the Fermi level, where the adsorption of MPcs leads to isolated molecules 16,5457 instead of the close-packed self-assembly observed on bulk Pb(111). This behavior has similarities with the present study where we found a reduction of the molecular packing on the Pb monoatomic layers. On coinage metals, the absence of dense assemblies was attributed to the residual charge accumulation at the molecule-substrate interface, resulting from an incomplete screening from the bulk electrons and leading to a subsequent electrostatic repulsion between molecules mediated by the screening of 2D surface states. Following this interpretation, we propose that our 2D Pb substrates may screen the residual molecular charges less eectively than the 3D Pb crystals thus promoting single and few-molecule clusters. Furthermore, on the SIC phase the presence of very dierent areas on the surface, i.e. anisotropic of the majority of molecules on specic zones of the phase giving rise to the observed frustrated molecular self-assembly. This mechanism might ultimately be further enhanced by the locality of the MnPc-substrate interaction leading to molecular adsorption only on certain surface sites and preventing dense molecular packing.Finally, it should be mentioned that MPcs have already been deposited on other 2D metallic systems. In particular, on the √ In phase of the In/Si(111) bilayer MPcs (M=Mn,Cu) self-assemble in closely packed square lattices. This indicates that, contrary to our observation, the molecule-molecule interaction is dominant in this system. Moreover, no signs of molecular-induced surface reconstruction were detected in the In/Si(111) bilayer. This fact highlights the peculiarity of our result and suggests that, as discussed below, in addition to the simple fact of having a low dimensional system, the intervention of the Pb/Si(111) interface might play an important role in tuning the balance between molecule-molecule and molecule-substrate interaction.

3 → 3 × √ 3 - 7 × √ 3 - 7 × √ 3 -

 3337373 SIC phase transition may ultimately be promoted by the energetic stability of the MnPc/SIC system, revealed by our DFT calculations, together with the high mobility of Pb atoms on top of the Si(111). This idea is reinforced by a recent study that reported about a charge transfer-induced structural phase transition of the √ Au/Si(111) monolayer under adsorption of 0.1 ML of Na atoms.[START_REF] Bondarenko | Large Spin Splitting of Metallic Surface-State Bands at Adsorbatemodied Gold/Silicon Surfaces[END_REF][START_REF] Bondarenko | Eect of Na Adsorption on the Structural and Electronic Properties of Si(111) √ 3 × √ 3-Au Surface[END_REF] ConclusionsIn conclusion, with this study we have shown that the adsorption of a minute amount of MnPc molecules (∼0.18 molecules/100 nm 2 ) on top of a Pb/Si(111) monolayer kept at 150 K and room temperature induces a structural transition of the √ Pb phase into the SIC-Pb phase on a macroscopic scale. Analysis of experimental and computational results revealed that the mechanism behind this surface transformation may be related to the strong charge transfer at the molecule/substrate interface as well as the strained nature of the surface reconstructions of the Pb monolayer. As shown by our DFT calculations, the phase transition could ultimately be promoted by the energetic stability of the MnPc/SIC-Pb/Si(111) system with respect to the MnPc/ √ Pb/Si(111) one.

√ 3 - 3 reconstructed 67 Acknowledgement 7 × √ 3 -

 336773 ter than 5×10 -11 mbar. Heavily n-doped Si(111) substrates were rstly outgassed at 873 K for several hours after being introduced in UHV and then repeatedly ashed at 1420 K by direct current heating. Atomically clean and defect-free Si(111)-7×7 reconstructed surfaces were prepared by repeatedly ashing the Si(111) sample to 1420 K and subsequent annealing from 1220 K to 820 K by slowly decreasing the temperature. Pb was evaporated at a rate of 0.65 ML/min on top of the freshly prepared Si(111)-7×7 surface kept at room temperature from an electron beam evaporator. A crystalline Pb monolayer was nally obtained by annealing the Pb/Si(111) system at about 600 K for a duration that depends on the desired surface reconstruction. In particular, after depositing 1.65 ML of Pb, the √ 7× Pb and the SIC phase were obtained by annealing the sample at ∼593 K for approximately 2 min 30 sec and 1 min 10 sec respectively. MnPcs were thermally sublimated at about 620 K on top of the Pb surface kept either at room temperature or at ∼150 K from a home-made resistively heated glass crucible.STM experimentsScanning Tunneling Microscopy (STM) measurements were performed in situ on a low temperature homemade system. Mechanically cut PtIr tips were used after being outgassed at around 470 K for several hours. All STM images were acquired in constant current mode with the sample kept at 77 K except that in gure 3 that was acquired at 300 mK. STM data were calibrated from high resolution images of the Si(111)-7×7 surface recorded at the same temperature.DFT calculationsTwo adsorption systems were considered in the calculations: a single MnPc molecule ad-Pb overlayer on Si(111) substrate respectively. Bulk Si(111) was modeled by means of a slab consisting of 6 Si bilayers passivated at one end with hydrogen atoms. In both adsorption systems, the atomic positions of the topmost Si bilayer, the top layer of Pb atoms and the MnPc molecule were all relaxed until atomic forces were smaller than 50 meV per Å. The projected spin-polarized density of states (PDOS) onto the d atomic orbitals of the Mn atom used a gaussian smearing with a width of 0.2 eV to smooth the plots. Spin-Polarized DFT calculations were performed for a periodic slab geometry as implemented in the VASP code. 64 A cut-o energy of 280 eV in the plane-wave expansion was used together with a 1 × 1 × 1 k-grid (Γ point) for the geometry optimization. For the static electronic structure calculations like PDOS and STM simulations a ner 3 × 3 × 1 k-grid was adopted together with an increased cut-o energy of 350 eV. The smearing technique using a Gaussian broadening of 0.01 eV was applied to achieve electronic convergence. Exchange and correlation were treated with the PBE form of the generalized gradient approximation. 65 Dispersion interactions were modeled with the variant of DFT-D3 method by Grimme. 66,We gratefully acknowledge F. Thibout (Kastler Brossel Laboratory) for the realization of the glass crucibles used during the experiments of this study. NL gratefully acknowledges nancial support from the spanish MICINN (project RTI2018-097895-B-C44) and FEDER funds. NW gratefully acknowlegde the nancial support of Labex Matisse and Investissement d'Avenir. CB gratefully acknowledges nancial support from the French ANR Rodesis (project ANR-16-CE30-0011-01). The computational work was granted access to the HPC resources of GENCI (TGCC and CINES) under grant n°A0050807364. Supporting Information Available Additional gures are available. They show the calculated atomic relaxation of the Pb monolayer after MnPc adsorption and the charge/spin density distributions of a MnPc molecule adsorbed on the √ Pb and the SIC-Pb phase.

  Pb phase, that covers the entire sample surface, is not a single domain phase but it consists of three equivalent domains 120 • mutually rotated and with an average size of about 100×100 nm 2 .

	transition.	2229 The	√	7 ×	√	3-

√

7 × √ 3-Pb surface before MnPcs deposition that develops below 270 K as a result of a reversible temperature-induced structural phase

  Figure 1b shows STM and LEED results of the Striped Incommensurate phase (SIC) which is the densest phase of the Pb monoatomic layer obtained at the saturation coverage of

	1.3 ML.	27,37 Here, it should be stressed the fact that, unlike the	√	7 ×	√	3-Pb phase, the SIC
	phase does not undergo structural changes from room temperature to low temperature.	27
	The analysis of STM images reveals that the incommensurate structure of this phase consists
	of small	√	3 ×	√	3-Pb domains (brighter regions) separated by meandering domain walls
	with quasi-	√	7 ×	√ 3-Pb atomic arrangements (dark stripes).

  Pb phase, our calculations reveal that the total molecular spin decreases from S=3/2 to S=1 as a result of charge transfer to the d xz /d yz orbitals that are hybridized with the ligand orbitals of the MnPc molecule.The S=1 state is ultimately given by the half-occupied d z 2 and d xy orbitals (see gure 5). It should be pointed out that the computed electronic structure of the MnPc adsorbed on the

	√ 7 ×	√	3-Pb monolayer is qualitatively similar to the one identied by Minamitani et al. 53
	for a MnPc adsorbed on Pb(111) or on a Pb island thick enough to mimic bulk properties.
	In contrast, the molecular magnetic moment is further reduced on the	√ 3 ×	√	3-Pb layer
	because the d z 2 orbital lls as a consequence of larger charge transfer. According to this the

d z 2 orbital should become non-magnetic. Accordingly, the total molecular spin is expected to be S=1/2 and arise from the only magnetic orbital that is left half-occupied, i.e. the d xy orbital. This orbital interacts weakly with the Pb/Si(111) monolayer and its occupation remains essentially unchanged, that is with roughly one unpaired electron.

In summary, our DFT analysis sheds light on the driving force of the MnPcs-induced

  Si atoms underneath the Pb layer (see gure S2 in Supporting Information). In this sense, the reduced dimensionality of the Pb monolayer permits a direct interaction of the MnPc molecule with the Pb/Si interface, leading to stronger and more localized bonds than on higher dimensional Pb substrates. For example, the recent calculations of Minamitani et al.

			15,43 In stark
	contrast to these results, our experiments reveal that MnPcs deposited on the Pb monolayer
	remain either isolated or self-assembled in trimers and their adsorption congurations reect
	the threefold symmetry of the SIC phase. This observation seems to suggest that the lower
	is the dimensionality of the Pb substrate the stronger becomes the molecule-substrate inter-
	action. This statement is supported by the values of adsorption energy (bonding distance)
	obtained by our DFT analysis, i.e. -2.86 eV (2.7 Å) and -2.74 eV (2.9 Å), which are higher
	(shorter) than those previously reported for a MnPc molecule adsorbed on bulk Pb(111),
	i.e. -1.52 eV	16 (3.29 Å	53

15,16,54

The small modulations of this weak molecule-substrate coupling arise from the incommensurability of a square lattice on the triangular Pb(111) surface leading MnPcs to lie on dierent adsorption sites. ). The computed spatial distribution of the electronic charge and spin density also underlines this strong interface eect: both densities accumulate heavily on the

  50,51 Jacob et al. 52 reported about the correlated electronic struc- ture of MnPc that retains its magnetic state, i.e. S=3/2, after adsorption on Pb(111) but that leads to an eective S=1 spin due to strong dynamical uctuations of charges between the d xz and d yz orbitals. Although our DFT results concerning the electronic/magnetic state of MnPc on the Pb monolayer are in qualitative agreement with the interpretation of Minamitani et al. Pb phase, in addition to the charge transfer to the d xz /d yz orbitals that reduces the total molecular spin by ∼1/2, there is an additional charge transfer that quenches the spin arising from the d z 2 orbital ultimately resulting in a total spin S=1/2 carried by the d xy orbital. This result further conrms the stronger interaction of MnPcs with the FePc or perylene. The explanation of this behavior was based on the fact that the herringbone surface reconstruction is itself induced by surface strain. According to this, the charge transfer that accompanies the molecular adsorption induces a change in the surface strain of Au(111) thus driving a new surface reconstruction. Therefore, adsorption of dierent molecules on Au(111) induces dierent surface patterns.

	MnPcs-induced surface reconstruction of the Pb monolayer
	√ While our work emphasizes a structural phase transition of the Pb monolayer from the 7× √ 3-
	53 for 3 to the SIC phase induced by MnPcs, no structural changes were ever reported on bulk Pb(111), those for the √ 3 × √ 7 × √ √ 3-Pb monolayer reveal something dierent. On the √ 3 × bulk Pb(111) or Pb/Si(111) islands. The reason of this may lie in the smaller strain at the √ Pb(111) surface compared to those within the Pb monolayer. To clarify this fact it is useful to mention that examples of surface modication induced by molecular adsorption can be 3-√ 3 × found on Au(111). Sun et al., 58 for instance, have reported modications of the herringbone √ 3-Pb phase than with the √ 7 × 3-Pb phase. √ surface reconstruction of Au(111) whose nal structure depends on the adsorbed molecule,
	i.e.
	Before concluding this section, we discuss the dierent charge transfers involving MnPc
	molecules on bulk or low-dimensional metallic substrates. On Au(111) the charge transfer
	proceeds from MnPc molecules to the substrate.	54 On Ag(111) it goes in the opposite di-
	rection	57 as found in our Pb monolayers and in 3D Pb(111). Despite this similarity, in the
	MnPc/Ag(111) system the charge transfer, i.e. 0.5 electrons,	57 is dominated by the organic
	ligands while on the	√	3 ×	√	3-Pb monolayer it is dominated by the Mn atom. This em-
	phasizes the locality of the interaction. Another interesting case is the adsorption of MPcs
	(M=Mn and Cu) on the 2D	√ 7× √ 3-In bilayer grown on Si(111). It revealed that the electron
	transfer goes from the substrate to the molecules in both cases.	19 However on one hand this
	charge transfer takes place on the organic ligands for CuPcs and yields 1.55 electrons.	19 On
	the other hand it is rather localized on the central atom in the case of MnPcs, in agreement
	with our ndings, and yields 1.78 electrons.	19
	The novelty of our result compared to the great variety of molecular behaviors encountered
	in literature until now is the new self-assembly of MnPcs on the Pb monolayer that re-
	ects a new way of tuning the molecule-substrate coupling by playing with the substrate
	dimensionality.

This material is available free of charge via the Internet at http://pubs.acs.org


Christophe Godin

Pascal Ferraro

Self-similar analysis of plant architecture reveals hierachical classes of meristem states

Keywords: Self-similarity, meristem, tree reduction, inverse problem

Introduction

Apical meristems are small embryogenic regions, located at the tip of plant axes, that build up plant organs by cellular division. The production of the meristems depends on their internal physical, physiological and genetic state and is controled by contextual factors (like micro-environment, availability of nutrients, etc.). In principle, the number of variables that may be used to define the state of a meristem, taking account the nature and the concentrations of molecules in each cell, their position, the physical stresses at each point, the geometry of cells, their genetic contents, etc., is infinitely large. Due to this intrinsic complexity, and to the current lack of hindsight on processes at such small scales, the connection between a meristem state, its micro-environment and what it produces at varying time scales seems until now largely out of reach.

However, the remarkable organization of plants at macroscopic scales makes the situation not so hopeless. The fact that plants are made up of the repetition of many similar components, at different scales e.g. [?, ?, ?, ?, ?, ?], provides macroscopic evidence for regularities and similarities in processes that drive meristem activity at microscopic scales.

In this paper, we propose to formalize this connection between macroscopic observations and microscopic, mostly invisible, processes. To achieve this connection, we formulate the following simplifying, though fundamental, scaling hypothesis : Scaling hypothesis: If two branching structures in a plant are similar, they were probably produced by meristems with similar state and context.

In this presentation, we shall show that it is possible to use this idealizedbut useful -hypothesis to organize the multitude of meristem states and contexts by classes of equivalence with respect to the similarity of what they produce.

In this first approach, the similarity between branching systems is considered to be purely structural (no geometry is taken into account for instance). Due to the nested nature of these structures, we show that the study of similarities between all parts of a plant boils down to studying the self-similar nature of the plant structure. Based on previous attempts to quantify self-similarity in plants [?, ?], we introduce a new method that enables us to define the degree of self-similarity of any plant.

This definition relies on the construction of a directed acyclic graph (DAG) corresponding to the reduction of the initial tree, when all the structural redundancy has been removed. Under the scaling hypothesis, the nodes of this graph can be interpreted as meristem states, and the edge between two states would denote the occurence of a meristem differentiation (from the source to the target state). Paths in the graph therefore denote all possible meristem differentiation sequences. If this DAG is linear (i.e. can be represented by a chain), the plant is said self-similar. In self-similar plants, there is thus a unique differentiation sequence for all the meristems of the plant. We show that it is possible to compute, for any given plant, the smallest self-similar plant which contains it. This makes it possible to define a unique differentiation sequence of meristem states for any plant and to associate each branching system of the original plant with one of these computed, hypothetical, meristem states.

The hierachical organization of the meristem computed in this way could be used as a guiding canvas to analyse the microscopic nature of meristems.

From biology to mathematical formalization: modelling plant architecture self-similarity Definition of an equivalence relation on branching systems. Different stratgies can be used to define an equivalence relationship between branching systems. They can be equivalent because they have the same root diameter, because they have the same size or because they bear the same number of flowers. Here we consider structural equivalence. Formally, this comes down to defining a notion of isomorphism between branching systems. In [?], we defined isomorphism between axial branching systems (i.e. branching systems for which a trunk is defined). Here we consider a less restrictive class of isomorphim between branching systems which holds for any type of tree structure.

Definition 1 (tree isomorphism) Let T 1 = (V 1 , E 1 ) and T 2 = (V 2 , E 2 ) be two rooted trees. A bijection φ from V 1 to V 2 is a tree isomorphism if for each (x, y) ∈ E 1 , (φ(x), φ(y)) ∈ E 2 .
To compute whether two branching structures T 1 and T 2 are isomorphic, we use a notion of edit-distance between trees [?, ?]. This distance is defined as the minimal number of elementary edit operations (insert, delete or match vertices) that is necessary to transform T 1 into T 2 : D(T 1, T 2) = min S si∈S s i . This distance has the following property:

Property 1 D(T 1, T 2) = 0 ⇔ T 1 ≡ T 1
Definition 2 (Reduction of a tree) Let T be a tree, we denote by R(T ) the graph obtained by quotienting T by the equivalence relation ≡.

It can be shown that this graph is a DAG and that there exists an algorithm that can compute this DAG in time O(|T | 2 ln|T |) [?]. A DAG for which there exists a path going through all its vertices is said linear.

Definition 3 (Self-Similar tree) A tree T is self-similar if R(T ) is linear. Among all the self-similar trees, we consider those trees that contain T , denoted as S(T ). Then, we consider trees T * in this set that minimize the distance with T .

Definition 4 (Smallest Self-similar Tree, SST) Let T be a tree and S(T ) be the set of all the self-similar trees that contain T . Then, we define the set of smallest self-similar trees containing T by: SST(T ) = {T * |T * = argmin T ∈S(T )

D(T, T )}

We show that, for any tree T , it is possible to compute an element T * ∈ SST(T ) in polynomial time O(), and give the corresponding algorithm.

Example: hierarchical organization of meristem states in Rice

The above approach was tested on different plant architectures. We present here results corresponding to the analysis of a rice panicle, already described in [?], Fig. 1. The topological structure of the panicle, T , is depicted in Fig 2 .  Figures ... We first compute the reduction tree S(T ), Fig. 3. This graph, from which the original tree can be reconstructed [?], is not linear and show a number of different meristem differentiation sequences. By computing an element T * of SST(T ) (Fig. 4), it is possible to find a single sequence of meristem state differentiation. The states of this sequence can be subsequently projected on the original topological structure, thus providing an interpretation of the entire structure in terms of meristem differentiation (Fig. 5)

Perspectives: from mathematical formalism back to biology

Due to the new technological and scientific progresses in molecular biology, imaging and developmental biology, the possibility to open up the meristem blackbox enables us to start seeking for processes and variables that determine at first order a meristem state and the nature of its growth context.