Ulisse Ferrari

Learning Maximum Entropy Models from finite size datasets: a fast Data-Driven algorithm allows sampling from the posterior distribution

Maximum entropy models provide the least constrained probability distributions that reproduce statistical properties of experimental datasets. In this work we characterize the learning dynamics that maximizes the log-likelihood in the case of large but finite datasets. We first show how the steepest descent dynamics is not optimal as it is slowed down by the inhomogeneous curvature of the model parameters space. We then provide a way for rectifying this space which relies only on dataset properties and does not require large computational efforts. We conclude by solving the long-time limit of the parameters dynamics including the randomness generated by the systematic use of Gibbs sampling. In this stochastic framework, rather than converging to a fixed point, the dynamics reaches a stationary distribution, which for the rectified dynamics reproduces the posterior distribution of the parameters.

We sum up all these insights in a "rectified" Data-Driven algorithm that is fast and by sampling from the parameters posterior avoids both under-and over-fitting along all the directions of the parameters space. Through the learning of pairwise Ising models from the recording of a large population of retina neurons, we show how our algorithm outperforms the steepest descent method.

Nowadays scientists from many different disciplines face the problem of understanding and characterizing the behavior of large multi-units complex systems with strong correlations [START_REF] Buzáki | Large-scale recording of neuronal ensembles[END_REF][START_REF] Phillips | A maximum entropy approach to species distribution modeling[END_REF][START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a population[END_REF][START_REF] Peyrache | Replay of rule-learning related neural patterns in the prefrontal cortex during sleep[END_REF][START_REF] Weigt | Identification of direct residue contacts in proteinprotein interaction by message passing[END_REF][START_REF] Cocco | From principal component to direct coupling analysis of coevolution in proteins: Low-eigenvalue modes are needed for structure prediction[END_REF][START_REF] Bialek | Social interactions dominate speed control in driving natural flocks toward criticality[END_REF][START_REF] Santolini | A General Pairwise Interaction Model Provides an Accurate Description of In Vivo Transcription Factor Binding Sites[END_REF]. Statistical Inference tackles these problems by inferring parameters of a chosen, context inspired, probability distributions to obtain models reproducing the system behavior. The basic strategy consists in choosing a model family described by a set of parameters and tune them to reproduce the dataset properties.

However, in many cases, due to a substantial unawareness of the system properties or to avoid biasing the results, hypotheses on the distribution functional form cannot be suggested nor trusted. To overcome this issue, the Maximum Entropy (MaxEnt) Principle [START_REF] Jaynes | On The Rationale of Maximum-Entropy Method[END_REF] suggests to search for the probability distribution with the largest entropy between those satisfying a set of constraints, which force the model distribution to reproduce the experimental averages of a list of observables. For systems of binary variables a common and fruitful choice is to constraint the model distribution to reproduce the experimental single and pairwise correlations. This choice has been successfully applied in system neuroscience [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a population[END_REF][START_REF] Cocco | Neuronal couplings between retinal ganglion cells inferred by efficient inverse statistical physics methods[END_REF][START_REF] Hamilton | Optogenetic Activation of an Inhibitory Network Enhances Feedforward Functional Connectivity in Auditory Cortex[END_REF][START_REF] Mora | Dynamical criticality in the collective activity of a population of retinal neurons[END_REF][START_REF] Tavoni | Inferred network of the prefrontal cortex activity unveils task-related coupling potentiations and cell assemblies[END_REF], gene regulation [START_REF] Lezon | Using the principle of entropy maximization to infer genetic interaction networks from gene expression patterns[END_REF], fitness estimation [START_REF] Ferguson | Translating HIV sequences into quantitative fitness landscapes predicts viral vulnerabilities for rational immunogen design[END_REF][START_REF] Mann | The Fitness Landscape of HIV-1 Gag: Advanced Modeling Approaches and Validation of Model Predictions by In Vitro Testing[END_REF] and many others fields of science. Moreover other possibilities with different observable lists have been also investigated [START_REF] Ganmor | Sparse loworder interaction network underlies a highly correlated and learnable neural population code[END_REF][START_REF] Tkacik | Searching for collective behaviour in a network of real neurons[END_REF], suggesting that depending on the context a careful choice of the observables can improve the inference accuracy and predicting power.

However, once the MaxEnt problem is posed as an inference task, finding its solution can be very hard. For large system size, in fact, the inference problem cannot be solved analytically and specifically devoted algorithms are required. The most known and widely used algorithm [START_REF] Ackley | A learning algorithm for boltzmann machines[END_REF] was introduced in the eighties and later developed [START_REF] Hinton | Training Products of Experts by Minimizing Contrastive Divergence[END_REF][START_REF] Broderick | Faster solutions to the inverse pairwise Ising problem[END_REF]. Other approaches include Selective Cluster Expansions [START_REF] Cocco | Adaptive cluster expansion for inferring Boltzmann machines with noisy data[END_REF][START_REF] Barton | Ising models for neural activity inferred via selective cluster expansion: structural and coding properties[END_REF], Minimum Probability Flow [START_REF] Sohl-Dickstein | New Method for Parameter Estimation in Probabilistic Models: Minimum Probability Flow[END_REF] and several approximation schemes [START_REF] Kappen | Efficient learning in boltzmann machines using linear response theory[END_REF][START_REF] Tanaka | Mean-field theory of Boltzmann machine learning[END_REF][START_REF] Aurell | Inverse Ising Inference Using All the Data[END_REF][START_REF] Ricci-Tersenghi | The Bethe approximation for solving the inverse Ising problem: a comparison with other inference methods[END_REF][START_REF] Jacquin | Efficient, fast and principled mean-field inference for strongly coupled data[END_REF].

In this paper we develop further the approach of [START_REF] Ackley | A learning algorithm for boltzmann machines[END_REF]. Our results are based on an analysis of the geometrical structure of the model parameters space. For the MaxEnt models inference this space is shown to benefit of peculiar properties, which allow us to introduce a novel quasi-Newton method, the Data-Driven algorithm, and to completely characterize its long-time learning dynamics. As it is affected by the randomness of Monta Carlo estimates, the parameters dynamics is stochastic and eventually converges to a stationary distribution around the log-likelihood maximum. The presented method takes advantage of this randomness and shapes the stationary distribution to reproduce the Bayesian posterior distribution and thus to sample from it. This last feature allows to avoid overfitting and endorses our algorithm to be well suited for dataset with highly inhomogeneous noise. We conclude with a test on biological data.

I. MAXIMUM ENTROPY MODELS

Systems of interest for MaxEnt approach are composed by N units that show a stochastic and coupled behavior. In order to be concrete, in this work we focus on binary, 0/1, units, but most of the results can be directly generalized to multiple-state variables, as Potts or Poisson models, or even continuous ones. In this framework, datasets are composed by B independent measurements of the synchronous state of the N system units: {{σ i (b)} N i=1 } B b=1 . As an example, for a binned spike trains recording of several neurons, σ i (b) = 1 may represent the activity (spike) of the i-th neuron in the b-th time-bin and σ i (b) = 0 its silence.

The first, crucial, step in the MaxEnt analysis is the choice of a set of observables. Observables are generic functions of the system units that should be chosen in order to catch the way system components interact. Strictly speaking, if some statistical feature is relevant for the system behavior, the corresponding observable should be included in order to force the model to reproduce that feature. On the other way round, if some feature is not essential, by excluding the corresponding observable, the model will adapt its behavior consistently with the other imposed features. Technically speaking the observables will be the sufficient statistics of the model probability distribution. As an example, the most common choice in the literature considers all the single variable terms (σ i ) and pairwise products (σ i σ j ) and it leads to the construction of the well know Disorder Pairwise Ising model [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a population[END_REF]. This particular choice allows to take into account the pairwise correlations between the units and allows the model to adapt higher order statistics consistently.

Moreover, a carefully choice of the observable should take into account the quality of the dataset. Noisy observables for which the average is not significantly estimated will induce data overfitting with the risk of strongly reducing the model prediction power. From this point of view and overloading the jargon of Bayesian Inference, the observables choice represents a sort of prior term (see sect. I B for details).

For the sake of generality, we like to present the results for an arbitrary observables choice and with this aim in mind we introduce a generic vector of observables Σ(σ) ≡ {Σ a (σ)} D a=1 , functions of the system units. The choice of the observables vector completely determines the functional form of the MaxEnt model family [START_REF] Jaynes | On The Rationale of Maximum-Entropy Method[END_REF][START_REF] Cocco | Adaptive cluster expansion for inferring Boltzmann machines with noisy data[END_REF]. In fact, by searching for the probability distribution which has the largest entropy among those that reproduce the observables averages, we obtain (see, for example [START_REF] Barton | Ising models for neural activity inferred via selective cluster expansion: structural and coding properties[END_REF], for the whole functional calculation):

P X (σ) = exp{Σ(σ) • X}/Z[X] (1)
where Z[X] is a normalization constant, X is a D dimensional fields vector, namely the model parameters, conjugated to the observable vector Σ and Σ(σ) • X = a Σ a (σ)X a is the scalar product in Euclidean space [30].

As the distribution (1) measures the probability of each possible system configurations, we can compute its dataset (log-)likelihood:

l data X ≡ B b=1 ln P X σ(b) . (2)
Even if later we will consider posterior sampling, see sect. I B, for the moment we restrict the inference task to the search for the set of fields X * that maximizes the loglikelihood:

X * ≡ arg max X l[X] (3)
where

l[X] = B X • P -ln Z[X] , (4)
and

P ≡ Σ ≡ 1 B B b=1 Σ(σ(b) ) (5)
are the experimental averages of the observables. In fact, finding the fields values X * solving

∇ a l[X * ] = B P a -Q a [X * ] = 0 , (6)
where ∇ a ≡ d/dX a , is equivalent to enforce the constraints:

Q[X * ] = P (7)
where

Q[X] ≡ Σ X ≡ Tr σ Σ(σ)P X (σ) . ( 8
)
are the model averages of the observables. Here and in the following . . . X means average over the model distribution (1) with fields X and (. . . ) always refers to dataset averages. Tr σ means summation over all possible system configurations.

For most of the reasonable observables choice and in case of good data quality, the solution of the maximization problem (3) exists and is unique. However for sake of completeness, in app. A we discuss when and how multiple solutions could arise.

A. The geometry of the X-space

In the following sections we will deal with the non trivial geometry embedding the fields space, the X-space. This geometry is described by three matrices: the (negative) log-likelihood Hessian H[X], the Fisher matrix I[X] and the model susceptibility matrix χ[X]

H ab [X] ≡ -∇ a ∇ b l X /B , (9)
I ab [X] ≡ ∇ a ln P X σ ∇ b ln P X σ X , (10)
χ ab [X] ≡ Σ a Σ b X -Σ a X Σ b X . (11)
H[X] describes the concavity of the log-likelihood function, I[X] describes the covariance of the log-likelihood gradient, whereas χ[X] that of the observables. If in a general optimization problem these three matrices differ, in the the MaxEnt model inference they coincide. In fact:

I ab [X] ≡ ∇ a ln P X σ ∇ b ln P X σ X (12) = -∇ a ∇ b ln P X σ X (13) = -∇ a ∇ b l X /B ≡ H ab [X] (14)
where in the first equality we use a well know properties of the Fisher matrix and in the second the fact that the log-likelihood is linear in P, the only data-dependent quantities. Moreover:

I ab [X] = -∇ a ∇ b ln P X σ X (15) = ∇ a ∇ b ln Z[X] (16) = Σ a Σ b X -Σ a X Σ b X ≡ χ ab [X] . (17) Indeed: χ ab [X] = I ab [X] = H ab [X] . (18)
The equations [START_REF] Tkacik | Searching for collective behaviour in a network of real neurons[END_REF] are the keystones of this study. They affect the geometry of the X-space and thus the inference in a peculiar way. In particular, the first equality allows us to introduce an efficient inference method, whereas the second allows us to completely characterize its long-time dynamics.

B. A Bayesian Framework for the Maximum Entropy Models

Until now, we posed the MaxEnt inference as a loglikelihood maximization problem, without considering that the finite size of the dataset could affect the estimate of the observables mean. The error in these estimates will inevitably result in some uncertainty on the fields inference that has to be taken into account. In fact, even if a carefully choice of the Σ avoids to include noisy observables, the system heterogeneity will results in a different precision in the fields estimation.

A Bayesian framework including prior and posterior distributions will exactly account for this uncertainty. In fact, through Bayesian inversion, we can compute the posterior distribution of the fields, P Post X data , from the prior distribution on the fields and the likelihood function, P data X = b P X σ(b) :

P Post X data = P data X P Prior X Norm . ( 19
)
where Norm is a X independent normalization constant. The width of the posterior distribution around its maximum quantifies the intrinsic uncertainty on the fields X and can be used to test the robustness of the inference. Explicitly, any scientific results based on a particular outcome of the fields inference should remain valid for all the fields sets with large posterior probability. Indeed the possibility to sample from the posterior is a powerful tool to test the robustness of the system analysis.

For the prior distribution we have two possible approaches: either include a flat distribution that does not depend on the fields value, either include a probability distribution that reflects some a priori knowledge on them. The first approach lets the dataset account for the whole uncertainly, and in case of very good dataset should be preferred. However, in applications dealing with strongly undersampled data and/or when some a priori knowledge is available, the second possibility has shown to be powerful.

In this work, as an example to show how to include a prior term, we focus on the L2-regularization of the form:

P Prior L2 X ≡ exp - B 2 X • η • X 2π |η|B D 2 , ( 20
)
where η is an arbitrary D-dimensional positive definite square matrix and |η| is its determinant.

II. VANILLA GRADIENT, NEWTON METHOD AND DATA-DRIVEN ALGORITHMS

Many of the algorithms suited for solving the MaxEnt inference task performs a dynamics in fields space that, starting from an initial condition, flows toward the maximum of the log-likelihood. In many applications, in fact, the log-likelihood gradient is fast to compute or estimate and can be used to drive the dynamics to the maximum. In this section we first review two of these approaches and then we propose the Data-Driven algorithm, the focus of this work.

A. Review: Vanilla Gradient and Newton Method algorithms

Before introducing the method proposed in this work, we like to review two well known inference method.

Ackley, Hinton and Sejnowski [START_REF] Ackley | A learning algorithm for boltzmann machines[END_REF] posed the inference problem as a dynamical process ascending the loglikelihood function along the gradient direction:

X t+1 -X t ≡ δX VG t = α(P -Q[X t ]), ( 21
)
where α is a learning rate. The algorithm works iteratively: at each time-step the computation of the model averages Q[X t ] allows to perform the fields update and to proceed towards the convergence, which is guaranteed for sufficiently small α. We call this approach Vanilla (Standard) Gradient (VG) algorithm.

Although it follows the gradient, VG will not go trough the shortest path even for arbitrary small α [START_REF] Amari | Natural Gradient Works Efficiently in Learning[END_REF]. The reason lies in the geometrical structure induced by the curvature of the log-likelihood function, namely its Hessian, see eq. ( 9). To take into account this geometrical effect we can multiply the log-likelihood gradient by the inverse of the Hessian [START_REF] Press | Numerical Recepies[END_REF] obtaining the well known Newton-Raphson method. However as suggested by Amari [START_REF] Amari | Why natural gradient?[END_REF][START_REF] Amari | Methods of information geometry[END_REF], in the curved manifold of the loglikelihood the natural local metric is the model Fisher matrix I[X], see eq. [START_REF] Cocco | Neuronal couplings between retinal ganglion cells inferred by efficient inverse statistical physics methods[END_REF]. Indeed, in the geometry induced by I[X], the steepest ascendant direction is I -1 ∇l, the contravariant form of the gradient (6) [START_REF] Amari | Natural Gradient Works Efficiently in Learning[END_REF]. However for the MaxEnt models inference, Hessian and Fisher matrix coincides and so do Newton-Raphson and Natural gradient methods. Consequently, we do not need to distinguish and simply replace the VG update [START_REF] Broderick | Faster solutions to the inverse pairwise Ising problem[END_REF] with

δX NM t = αχ -1 [X t ] • (P -Q[X t ]) (22)
to obtain both the Newton Method (NM) and the Natural gradient. The positiveness of χ ensures the convergence of this method at least for infinitesimally small α (see later). Despite it is optimal, the Newton Method is slowed down by the time-consuming estimation and inversion of χ[X] at each update step. In the following we suggest a way to bypass this problem.

B. The Data-Driven Algorithm.

As χ[X] depends only on model averages of observables products, we can approximate its value at the solution X = X * with the dataset configurations list:

χ ab [X * ] ≈ χ ab ≡ Σ a Σ b -Σ a Σ b . (23)
This matrix can be computed before running the inference dynamics and then used to evaluate the fields update. The resulting Data-Driven (DD) quasi-Newton Method update rule reads:

δX DD t = α χ -1 • (P -Q[X t ]) . ( 24
)
The quality of the approximation χ[X * ] ≈ χ depends on two main hypotheses: i) the ability of the MaxEnt model to reproduce dataset statistical properties beyond the mean of Σ and in particular the experimental Σcovariance and ii) the good sampling quality of the experimental dataset. The first hypothesis, however, was already partially assumed when it has been chosen the observables to reproduce and so the MaxEnt model to fit. Otherwise stated, if the approximation is poor, it means that the chosen MaxEnt model was not a good choice. The second hypothesis, instead, reflects the quality of the whole inference task. In case of strong undersampling the inference problem, despite being mathematically well posed, is meaningless as the information encoded in the dataset does not support the fine tuning of the fields. In sect. VI we will give a practical condition to test this second hypothesis.

The quality of the approximation (23) controls the speed up factors of the DD algorithm: the better the approximation is, the better and faster the DD will work. On the contrary, when the approximation is poor, the DD will not be advantageous with respect to the VG. In fact, as χ is positive definite, no matter how bad the approximation is, for sufficiently small α the DD will still converge toward the right solution X * .

In conclusion, for cases where the inference problem is meaningful and the consequently the approximation ( 23) is valid, the DD algorithm will speed up the inference, whereas in the other cases DD will be useless but not counterproductive.

III. THE LEARNING DYNAMICS

All the approaches explained before, (VG, NM and DD) solve the learning task through a discrete-time dynamics in the fields space. At each time step the estimation of the model averages of the observables (Q[X]) allows to compute the fields update and continue the dynamics. For large system size N , however, the exponential complexity of the problem prevents the exact computation of the observables averages and some approximations are required. A standard, but fruitful, choice consists in using Markov-Chain Monte-Carlo (MC) with Metropolis algorithm to sample M system configurations and use them to estimate:

Q MC X ≡ 1 M M b=1 Σ(σ MC (b) ) . ( 25
)
Q MC X is now a random variable approximating Q[X] up to O( 1/M ) fluctuations.
Even for large M , the randomness of Q MC X will eventually affect the convergence of the dynamics, as for X t sufficiently close to X * , the size of the gradient will become comparable with its fluctuations. However, at the beginning, we expect the dynamics to be almost deterministic and then to become stochastic only at the end. We indeed separate the dynamics in two regimes:

1. Approaching the convergence, when P -Q MC X is small but still much larger than its fluctuations.

2.

The long-time stochastic dynamics, when P -Q MC X 1/M and thus comparable with its fluctuations.

In appendices B and C we will provide several details of the two regimes and here we only present the mayor results.

A. Approaching the convergence

After an initial transient where the dynamics strongly depends on the chosen algorithm (VG,NM or DD) and on the initial conditions, we expect the X t to approach the log-likelihood maximum X * , so that we can approximate the log-likelihood function up to the quadratic order. Given δl

[X] ≡ l[X] -l[X * ], we have δl[X] - B 2 (X -X * ) • χ[X * ] • (X -X * ) ≈ - B 2 (X -X * ) • χ • (X -X * ) . ( 26
)
In this approximation the dynamics is exactly solvable upon projecting the fields on the χ[X * ] Eigenvectors {V µ } D µ=1 : δX µ t = a V µ a δX a,t . Along a µ-Eigenspace the convergence of the VG algorithm scales with the corresponding Eigenvalue λ µ as:

δX µ t ∼ (1 -αλ µ ) t .
Consequently, to ensure the algorithm convergence, we need α < 2/λ µ along all directions. Moreover, by optimizing the convergence speed along all directions simultaneously we obtain α

(VG) BEST = 2/(λ > + λ < )
, where λ >/< are the largest/smallest Eigenvalue. In App. B we present some details and here we simply notice how α

(VG)
BEST can be squeezed to very small values by a large λ > preventing the learning along all the direction with λ µ λ > . Consequently, for dataset where the Eigenvalues of χ[X * ], or of its approximation χ , spread over several order of magnitude the convergence of the VG algorithm will be very slow.

Because of the quadratic approximation, the NM and the DD algorithms coincide and we do not distinguish between them. Within their dynamics, δX µ t ∼ (1 -α) t independently of the Eigenvalues λ µ . consequently α < 2 is the only convergence condition and α (DD) BEST = 1.

B. The long-time stochastic dynamics

Before proceeding we like to introduce the shortcut notation:

N [m; s](t) ≡ e -1 2 ab (ta-ma)s -1 ab (t b -m b ) (2π) D |s| , (27)
a normal distribution with average m and covariance s evaluated at t. When X → X * , we expect ∇l MC X ≡ B P -Q MC X → 0 only on average, with not negligible fluctuations of O( 1/M ) . In this regime the fields dynamics is not anymore deterministic and the convergence becomes a stochastic process. Consequently, rather than converge to a fixed point, the fields will approach an equilibrium stationary regime around X * .

The stochastic process is ruled by the discrete-time master equation:

P t+1 (X ) = DX P t (X) W X→X , (28)
where the transition rates depend on the distribution of ∇l MC X .

In the large M limit, the gradient distribution can be approximated as a Gaussian to obtain an analytic expression of W X→X , see app. C. By asking P t (X) to be invariant under the evolution [START_REF] Ricci-Tersenghi | The Bethe approximation for solving the inverse Ising problem: a comparison with other inference methods[END_REF] we can obtain the stationary distribution P ∞ (X):

P VG ∞ (X) = N X * ; α M 2δ D -α χ -1 (X), (29)
P DD ∞ (X) = N X * ; α M (2 -α) χ -1 (X), (30)
where δ D is the identity matrix in dimension D. Here the typical fluctuations of X around X * must consistently verify the approximation X ≈ X * : (X -X * ) 2 P∞ should be small enough to allow the expansion [START_REF] Tanaka | Mean-field theory of Boltzmann machine learning[END_REF].

In the stationary regime, on top of the fluctuations induced by the MC, the X distribution will induce a second source of noise in the actual distribution of Q MC . On average we expect (see app. C):

P VG ∞ (Q MC ) = N P; 2 χ (2δ D -α χ ) -1 M (Q MC )(31) P DD ∞ (Q MC ) = N P; 2 χ M (2 -α) (Q MC ) . (32)
Interestingly, the two algorithms provide different Q MC distributions. In particular, the VG fluctuations along directions with large (small) λ µ are larger (smaller) than the DD ones. These differences will have consequences on the ability of both algorithms to reproduce the dataset statistics and thus to avoid both under-and over-fitting.

C. The dynamics under an external stochastic force

For forthcoming purposes, here we characterize the learning dynamics when a linear stochastic force term is added to the gradient term in the learning rules. This analysis will be useful when a L2 prior term is included in the inference procedure. We modify eq. ( 24) as:

δX DDη t = α χ η -1 • ∇l MC X + F η X . ( 33
)
where

χ η ≡ χ + η (34)
P F η X = N -η • X; η M (F η X ). ( 35
)
The calculation for the stationary fields distribution follows as in the previous section and it results in:

P DDη ∞ (X) = N X * η ; α M (2 -α) χ η -1 (X), (36)
where

X * η = χ η -1 • χ • X * . (37)
Analogously to [START_REF] Press | Numerical Recepies[END_REF], from (36) it follows:

P DDη ∞ (Q MC -F η ) = N P; 2 χ η M (2 -α) (Q MC -F η ) .
(38) Expressions for the VG algorithm can be easily obtained from [START_REF] Machta | Parameter Space Compression Underlies Emergent Theories and Predictive Models[END_REF] and [START_REF] Peyrache | Spatiotemporal dynamics of neocortical excitation and inhibition during human sleep[END_REF] with substitution α → α χ η -1 .

IV. AVOIDING UNDER-AND OVER-FITTING BY SAMPLING FROM THE POSTERIOR DISTRIBUTION.

In the previous section we shown how any actual implementation of the algorithm is stochastic. In particular, depending on the chosen algorithm and its parameters, α and M , we expect a whole probability distribution on the output fields X. This stochasticity raises questions on which implementation should be preferred in practical applications. The simple strategy of trying to obtain the best possible approximation of X * by increasing M or reducing α will unambiguously lead to data over-fitting thus limiting the model prediction power. The presence of noise in the estimation of P induces some uncertainty on the model fields X that should be taken into account to avoid over-fitting.

This fields uncertainty is quantified by the posterior distribution of the fields given the data, see eq. [START_REF] Ackley | A learning algorithm for boltzmann machines[END_REF].

Within the approximation [START_REF] Tanaka | Mean-field theory of Boltzmann machine learning[END_REF], the posterior simplifies to

P Post X data ∝ e -B 2 (X-X * )• χ ab •(X-X * ) (39)
in the case of flat prior and to

P Post X data ∝ e -B 2 [ (X-X * )• χ •(X-X * )+X•η•X ] ∝ e -B 2 (X-X * η )• χη•(X-X * η ) (40)
when an L2 prior, see eq [START_REF] Hinton | Training Products of Experts by Minimizing Contrastive Divergence[END_REF], is considered (X η and χ η are those defined in eqs. ( 37) and ( 34)).

If by tuning the algorithm implementation and settings we can match the a stationary fields distribution and the posterior, we will be able to sample from it, thus avoiding any over-and/or under-fitting. In case of flat prior, we have to compare [START_REF] Welling | Bayesian learning via stochastic gradient langevin dynamics[END_REF] with ( 29) and (30), whereas in case of L2 prior (40) with [START_REF] Machta | Parameter Space Compression Underlies Emergent Theories and Predictive Models[END_REF] for the DD algorithm and the analogous for the VG one: by setting

α = 2M B + M (41)
the DD algorithm will sample from the correct posterior distribution and over-and under-fitting will be avoided along all the D dimensions. On the contrary, for the VG algorithm such a setting does not exist and the algorithm will not sample from the posterior. In particular directions with αλ µ 2M/(B + M ) are strongly over-fitted and those with αλ µ 2M/(B + M ) are strongly underfitted. Otherwise stated: by choosing the VG algorithm with α small enough to well reproduce directions with large λ µ will result in overfitting of all directions with small λ µ .

Among the consequences of overfitting, for finite B we expect an overestimation of the observed model loglikelihood. In fact, as the exact inference fields X * reproduce also the noise in the experimental averages P, we expect l * ≡ l[X * ] > l, where l is the log-likelihood of the model that generates the data, the true log-likelihood. In order to quantify this effect, in appendix D we compute the average log-likelihood estimation in the case of the exact inference when the data are synthetically generated by a MaxEnt model with true fields X. By averaging over the distribution of P, we find that

l * -l = D 2B , (42)
which shows how the log-likelihood maximization, see eq.

(3), induces a finite bias leading to a log-likelihood overestimation. In the case of the DD algorithm, instead, the average over the posterior distribution exactly cancels the bias, and the true log-likelihood value is recovered.

V. ALGORITHM IMPLEMENTATION

The core of DD algorithm is to iteratively update the fields X t with the rule [START_REF] Sohl-Dickstein | New Method for Parameter Estimation in Probabilistic Models: Minimum Probability Flow[END_REF], where Q[X t ] is approximated by Q MC X through a MC sampling of M configurations. Ideally α = α DD BEST = 1 and M = B is the fastest setting that satisfies the condition (41). However any practical implementation will face two main difficulties. The first lies in the fact that the quadratic approximation [START_REF] Tanaka | Mean-field theory of Boltzmann machine learning[END_REF] may not be valid. This will happens at the beginning of the dynamics, when X t if far from X * , but also when B is not large enough to have (X t -X * ) 2 small enough to discard third order terms. If it is the case, the distribution (30) will be non-stationary. For this reason we allow the algorithm to adapt the value of α at each time-step. The second difficulty reflects the fact that the algorithm is not suppose to converge to X * but rather to a probability distribution. Consequently we need a condition that signals the onset of the thermalization and allows us to start storing X t s as samples from the posterior distribution. For this reason we introduce the following quantity:

t ≡ B 2D P -Q MC Xt • χ -1 • P -Q MC Xt . (43)
Under the distribution (32) with α = 1 and M = B, ∞ will have distribution:

P (ε ∞ ) = 2 D D 2 Γ( D 2 ) 2 D 2 ε D-1 ∞ exp - D 2 ε 2 ∞ ( 44
)
where Γ(x) is the Gamma function. ∞ has:

ε ∞ = 2 D Γ( D+1 2 ) Γ( D 2 ) ----→ D→∞ 1 (45) ε 2 ∞ -ε ∞ 2 = 1 -ε ∞ 2 ----→ D→∞ 0 . ( 46
)
In the large D limit, the ∞ -distribution shrinks to a Dirac-delta function at ∞ = 1: if before thermalization we expect with high probability ε t 1, once the algorithm starts sampling from the stationary distribution (30) we expect ε t ≈ 1. Consequently, once the condition ε t ≤ 1 is full filled, the subsequent X t are good estimations of the fields. In order to effectively sample from (30) it will be still necessary to keep running the algorithm in order to decorrelate from the initial condition X 0 .

The DD algorithm that we implemented can be sketched as follows:

1. Chose an initial configuration for X 0 compute/evaluate Q[X 0 ], see eq. ( 8) and compute ε 0 from eq. ( 43). Then set α 0 = 1 and M 0 = min( B ε 2 0 , B).

Iterate the following step:

(a) update the X t through eq. ( 24),

(b) estimate Q[X t ] though M t = min( B ε 2 t-1 , B)
Gibbs samplings, see eq. ( 25), (c) compute ε t from eq. ( 43) and i. if ε t < ε t-1 accept the update and set α t = α t-1 δ + ii. otherwise discard the update, set α t = α t-1 /δ -and estimate again Q[X t ].

3. As soon as the condition t < 1 is full filled, (a) either fix α and let the distribution P t (X) decorrelate from X 0 and thermalize to P ∞ (X), see ( 30). (b) either stop the algorithm and retain X t as a fields list solving the inference problem.

A variable and adapting α is required because the system can be be far outside the validity range of [START_REF] Tanaka | Mean-field theory of Boltzmann machine learning[END_REF]. In order to avoid cycles, we heuristically set δ + = 1.05 and

δ -= √ 2.
The choice of an adapting M t = min( B ε 2 t , B) ≤ B, instead, allows to save time during the deterministic dynamics regime, see sect. III A, when the algorithm does not require a high precision in the gradient estimate. ε t , in fact measures also the norm of the loglikelihood gradient in the metric defined by its fluctuations,

B 2 I[X]/M = B 2 χ[X]/M ≈ B 2 χ /M , see eq. (C4): ε t = M B ∇l[X t ] B 2 χ /M . (47)
For M t = B/ε 2 t we expect the estimate of the gradient to be statistically confident within one standard deviation and consequently enough well estimated. Moreover the randomness induced by small M t values decorrelates the dynamics from the initial condition. When the algorithm approaches convergence, ε t → 1 and consequently M t → B thus full filling condition (41).

The complexity of the algorithm depends linearly on N B through the MC sampling. From the tests we performed the number of required steps depends mostly on the dataset properties and in particular on the exactness of χ[X * ] ≈ χ . The hardest limit of the DD lies in the memory allocation for storing χ . Working in double precision, for an hardware with 32Gb of RAM the maximum number of manageable units is N ∼ 350 [START_REF]As the matrix χ has several symmetries a compressed encoding could[END_REF].

A. The L2 prior term

In order to include a L2 prior term we add to the dynamics the stochastic force term introduced in sect. III C. At each iteration, we should modify the update rule by adding to the gradient term ∇l MC X an independent realization of the random variable F η X according to the distribution [START_REF]As the matrix χ has several symmetries a compressed encoding could[END_REF]. As the prior term will biases the inference we need to modify the function (43) in order to account for the difference between ( 32) and ( 38): 23), for a pairwise Ising model inference problem on biological data. As often happens in biological problems [START_REF] Machta | Parameter Space Compression Underlies Emergent Theories and Predictive Models[END_REF], the Eigenvalues span several order of magnitude enhancing the heterogeneity of the parameter space. As will be show in the inset of Fig. 3, the susceptibility matrix of the inferred Ising model has a very similar spectrum. Consequently, even with the optimal learning rate, α

η t = B 2D P -Q + F η Xt • χ η -1 • P -Q + F η Xt (48)
(VG)
BEST of eq. ( B3), VG will be very slow in learning along the directions with the smallest Eigenvalues. As the sensitivity (1/B) lies well below the smallest Eigenvalue all the susceptibility spectrum is informative and has to be reproduced by the inference. See text for more details. Data from a ∼ 2.1h (binned at 16ms) recording of 95 rat retinal ganglion cells subject to visual stimulation [START_REF] Mora | Dynamical criticality in the collective activity of a population of retinal neurons[END_REF].

which, remarkably, does not require the explicit inversion of the matrix χ . Moreover, as χ is, by construction, non-negative, χ η is a positive matrix and can be inverted. As a consequence, the DD algorithm regularized with L2 prior can be applied when an undersampling induces zero modes in the empirical estimation of the susceptibility matrix.

VI. TEST

As explained before, the DD algorithm performance depends mostly on how good is the chosen MaxEnt model in modeling the dataset statistical properties. As expected, we succeed in inferring back the fields X from synthetic dataset obtained through simulation of Max-Ent models. Indeed we find more interesting to present an application to biological data where the underlying statistical distribution does not belong to class of models considered in the Inference task.

We tested the DD algorithm on an ex-vivo multielectrode array recording [START_REF] Marre | Recording of a large and complete population in the retina[END_REF] of 95 rat retinal ganglion cells [START_REF] Mora | Dynamical criticality in the collective activity of a population of retinal neurons[END_REF]. The retina was stimulated though a video showing two randomly moving bars [START_REF] Mora | Dynamical criticality in the collective activity of a population of retinal neurons[END_REF] displayed with a frame rate of 60Hz. The ∼ 2.1h spike trains recording was binned at 16ms to obtain B ∼ 4.8 • configurations, where we assign σ i (b) = 1 if cell i emitted at least one action potential in the time-bin b and σ i (b) = 0 otherwise. We infer a pairwise Ising model thus restricting the observables list to single and pairwise correlations ({Σ a } D a=1 = {{σ i } N i=1 , {σ i σ j } N i<j=1 }) and consistently the fields to biases and pairwise interactions

({X a } D a=1 = {{h i } N i=1 , {J ij } N i<j=1 }).
With this observables choice the model takes the form of the well known J DDη -0.5 0 0.5 Blue triangles compare the inferred fields with (x-axis) and without (y-axis) prior term. Red circles compare the inferred fields with regularization (x-axis) with the expected ones (yaxis):

X DD η ≡ χη -1 • χ • X DD , see eq. ( 37
). The black line shows equality. As a consequence of the L2 regularization the inferred fields are smaller in absolute value. Dataset as in fig. 1.

Disordered Ising Model:

P h,J (σ) = exp i h i σ i + i<j J ij σ i σ j / Z[h, J] (49)
where Z[h, J] is the normalization constant.

Before performing the inference it is important to check whether the the number of measurements is large enough to empirically estimate the observables averages with good precision. In the large B limit, we expect an empirical error δP in the estimate of P to have zero mean and covariance χ /B, see eq. (D2). By linear regression, en error δP in the data will induce en error:

δX = δQ[X] δX X=X * -1 • δP ≈ χ -1 • δP (50)
on the inferred fields. Consequently, by approximating χ ≈ χ :

P (δX) = N 0 ; (B χ ) -1 ](δX) . (51)
To ask for this error to be small, we should indeed require:

B λ µ 1 , for all µ ⇒ 1 B λ MIN . (52)
As shown in Fig. 1, in the considered dataset, condition ( 52) is satisfied: the experimental sensitivity, 1/B is much smaller then the smaller Eigenvalue of χ . On an Intel Core TM i7-4770 with eight cores at 3.4GHz the written in Matlab DD algorithm takes 100 ± 5s to reach convergence. In Fig. 2 we show the behavior of both t and α t against the number of inference steps. As the number M of MC samplings is not constant during the inference, steps are of different time duration.

In order to show the quality of the inferred fields, in the main panel of Fig. 3 we scatter-plot the experimental connected correlations against those estimated through B sampling of the inferred Ising model distribution. In order to give an insight on the validity of [START_REF] Barton | Ising models for neural activity inferred via selective cluster expansion: structural and coding properties[END_REF], in the inset we scatter-plot the ordered Eigenvalues of χ (xaxis) against those of χ[X * ].

To convince the reader that the tested inference problem was difficult, we apply two example Mean-Field approaches to the dataset. In Fig. 4 we show the scatterplot of the inferred fields against those of Naive [START_REF] Kappen | Efficient learning in boltzmann machines using linear response theory[END_REF] and Resummed [START_REF] Jacquin | Efficient, fast and principled mean-field inference for strongly coupled data[END_REF] mean-field approximations. Despite Resummed works much better than Naive, for this dataset both approximations largely overestimate couplings and biases.

For comparison, we tested the VG algorithm on the same dataset with M = B and various α ≤ α (VG)

BEST . With α = 0.2α (VG)
BEST the VG takes ∼ 4.2 • 10 4 s to converge, whereas for larger α values it was not able to satisfies the early-stopping condition < 1.

In Fig. 5 we scatter plot in blue the inferred fields against those inferred with a L2 regularization term, see eq. ( 20), with η = 5 • 10 -3 δ D , whereas in red we scatter the inferred with prior fields against the expected ones, see eq. [START_REF] Marre | Recording of a large and complete population in the retina[END_REF].

A part from the retina example, we succeeded in inferring the Ising model fields from synthetic data, human temporal cortex recording [START_REF] Peyrache | Spatiotemporal dynamics of neocortical excitation and inhibition during human sleep[END_REF] , rat pre-frontal cortex recording [START_REF] Peyrache | Replay of rule-learning related neural patterns in the prefrontal cortex during sleep[END_REF] and others. In all the tested cases convergence times are of the order of tens to few hundreds of seconds.

VII. CONCLUSIONS AND DISCUSSION

In this study we introduced a Markov-Chain Monte-Carlo based algorithm for solving the MaxEnt inference problem and sampling from the posterior probability distribution of the MaxEnt fields X, namely couplings and biases of an Ising model or their generalization in the presence of non-pairwise interactions. We carefully analyze the learning dynamics and separate two different regimes: i) a deterministic dynamics that approaches the solution and ii) a stochastic dynamics of the probability distribution of the fields, P t (X), that thermalizes to a stationary distribution. By tuning the algorithm settings, namely α and M , this distribution can reproduce the posterior distribution of the inference fields thus allowing to sample. We concluded presenting an implementation of the algorithm and a test on a biological dataset showing how the presented algorithm outperforms the standard learning approach. The core of the algorithm is the approximation ( 23) which requires at first to have enough data to properly estimate the empirical susceptibility then that the probability model chosen for the inference reproduces quite well the data statistics. As the applications to biological data have shown, both conditions have not to be intended as rigid constraints, but rather as requirements for obtaining the largest advantage of the DD approach.

The key properties of the MaxEnt inference that driven this study are the relationships [START_REF] Tkacik | Searching for collective behaviour in a network of real neurons[END_REF]. The equality among the log-likelihood Hessian and the model Fisher allows to tune the MC fluctuations of the log-likelihood gradient (ruled by the Fisher) in order to reproduce the fluctuations of Posterior distribution around the log-likelihood maximum (ruled by the Hessian). The equality with the susceptibility matrix allows to obtain an expression of both the Fisher and Hessian that depends only on Xdependent averages of X-independent functions that at the solution X * can be approximated through the data.

In sect. III, we have carefully analyzed the learning dynamics and in III B we have characterized it as a stochastic process that converges to a stationary distribution. Later, in sect. IV, we have shown how the DD algorithm stationary distribution can be tuned to the posterior distribution. This characterization has several advantages:

• the computation of P DD ∞ (Q MC ), together with the running evaluation of t , see eq. ( 43) provides an useful early-stopping condition for the algorithm that allows to check when the inference has been accomplished. Consequently it is useful to save much computational time.

• any inference algorithm has to account for the randomness of the MCMC fluctuation and typically large computational efforts are required to get rid of this noise source (large number of MCMC samplings M ) . The DD algorithm, instead, takes advantage of these fluctuations to avoid overfitting and to decorrelate from the initial condition.

• as we have shown in appendix D, in comparison with the exact inference, the posterior sampling allows to avoid overfitting the model.

• as in the case of the VG algorithm, an overfitting along several directions happens with all the algorithms that do not rectify the fields space.

• as explained in sect. I B, the posterior sampling can provide several consistent lists of inferred fields, thus allowing to test the robustness against the dataset noise of any analysis based on the inferred fields.

Our approach and the use of the MC randomness to induce noise in the algorithm outcome may remind the application of stochastic gradient method for sampling from the fields posterior distribution [START_REF] Welling | Bayesian learning via stochastic gradient langevin dynamics[END_REF]. Despite the results may look similar, in order the induce fluctuations, here we take advantage from the intrinsic algorithm randomness instead of introducing it by sub-sampling the dataset. The mayor advantage results from the possibility to select and optimize properly the algorithm setting in order to get rid of the MC induced noise that otherwise will affect the dynamics in a spurious way. One important drawback arises from the fact that the algorithm fluctuations may not reproduce the experimental ones when the Gaussian approximation ( 26) is not valid and second order terms are not enough to approximate the log-likelihood function around the solution. In this case, however, one can notice how the equalities ( 18) can be extended up to higher order cumulants. These relationships, together with an appropriate update rule that takes into account higher order corrections may extend the equivalence between P ∞ (X) and the posterior. We let this generalization for forthcoming investigations.

The DD approach takes its place in the list of algorithms for exactly solving the MaxEnt model inference problem. Its greatest strengths are the velocity and the possibility to sample from the posterior, whereas its stronger limitation is the memory requirement for the storage of χ . As an example, the inference of a N = 350 pairwise Ising model requires 32Gb of RAM.

Depending on the dataset properties the DD should or should not be preferred to others algorithms as [START_REF] Broderick | Faster solutions to the inverse pairwise Ising problem[END_REF], [START_REF] Cocco | Adaptive cluster expansion for inferring Boltzmann machines with noisy data[END_REF] or [START_REF] Sohl-Dickstein | New Method for Parameter Estimation in Probabilistic Models: Minimum Probability Flow[END_REF]. In general, if the condition (52) is largely satisfied, we expect DD to be the best choice, because the approximation ( 23) is expected to be valid. In the opposite case, a largely unsatisfied (52) suggests that the inference of the chosen MaxEnt model is not meaningful and the observables list has to be modified. However for cases in between some tests with different algorithms have to be performed. In particular cases, some available a priori knowledge of the system could help in the algorithm choice. For example, if the underlying interaction graph is naturally clusterized in almost non-interacting subcomponents, Selective Cluster Expansion (SCE) [START_REF] Cocco | Adaptive cluster expansion for inferring Boltzmann machines with noisy data[END_REF] is probably the best choice. SCE, in fact, splits the system in many subunits that are recursively joined together to form larger and larger building blocks of the reconstructed interaction network. If the interaction graph is clusterizable, SCE will recognizes these units and accomplish the inference task quickly. However, up to our knowledge, a simple and generic argument for choosing the best suited algorithm for the actual inference problem is still missing and it would be of large interest.

of very low convergence. The learning rate that optimize the convergence speed along all direction simultaneously is and the solution reads:

α (VG) BEST = 2 λ > + λ < , (B3)
In particular, α

(VG)
BEST can be squeezed to small value by large λ > preventing the learning along the direction with λ µ λ > . As an example, for the biological data we will consider, see Fig. 1, λ > /λ < 10 5 . Moreover the ratio λ > /λ < has been shown to diverge in synthetic data of model at criticality [START_REF] Machta | Parameter Space Compression Underlies Emergent Theories and Predictive Models[END_REF].

λ MAX or α 2/λ MAX it can reach very large values, thus nullifying the inference effort.

Note that the three matrix appearing in the covariance of eq. (C1) and in the mean and covariance of eq. (C4) are a priori different: the first is the model susceptibility, the second is the log-likelihood Hessian where the third is the model Fisher matrix. However, as explained in section I A for the MaxEnt inference problem these three matrices coincide providing the results ( 29) and (30). J ij = 0, , otherwise h i = -0.46 + j =i J ij (D8)

We first computed l = 5.7092 and then we simulated the model S = 2 • 10 4 times to collect S empirical estimates of δP. The number of MCMC sampled was fixed at B = 2 13 . For each of these S realization we estimate l * δP and l DD δP and in Fig. 6 we compare their histograms. Dotted lines correspond to the theoretical distributions (D4) and (D7), whereas the arrows indicate the empirical means. The black vertical lines represents l. As can be appreciated by the small difference between the black vertical bar and the blue arrow, the average over the posterior distribution distribution removes the bias produced by the exact inference. Moreover the difference between the red and blue arrows is approximately D/(2B), equals to the length of the horizontal segment. More precisely:

FIG. 1 .

 1 FIG. 1. (Colors online) Eigenvalues distribution of the empirical susceptibility matrix, see eq. (23), for a pairwise Ising model inference problem on biological data. As often happens in biological problems[START_REF] Machta | Parameter Space Compression Underlies Emergent Theories and Predictive Models[END_REF], the Eigenvalues span several order of magnitude enhancing the heterogeneity of the parameter space. As will be show in the inset of Fig.3, the susceptibility matrix of the inferred Ising model has a very similar spectrum. Consequently, even with the optimal learning rate, α

FIG. 2 .FIG. 3 .

 23 FIG. 2. (Colors online)Behavior of the pairwise Ising model inference of biological data through the Data-Driven algorithm. εt (lower blue curve, left axis) and of αt (upper red curve,right axis) during the algorithm running plotted against the iterative step number. Abscissa is not proportional to the running time as each step performs a different number M of MC sampling, see text. αt is raised (lowered) at every step where εt decreases (increases). The algorithm requires 100 ± 5s to produce a set of fields X satisfying the early-stop condition εt < 1. Dataset as in fig.1.

FIG. 4 .

 4 FIG. 4. (Colors online) Scatterplots of the inferred couplings (main panel) and biases (inset) with the DD algorithm against those of Naive (blue triangles) and Resummed (red circles) Mean-Field. The black lines show equality. All the approximations fails in the reconstruction with a tendency to overestimate both couplings and fields (notice the different axis scale). Dataset as in fig. 1.

  FIG. 5. (Colors online) Scatterplot of the inferred couplings (main panel) and biases (inset) for a pairwise Ising model with and without L2 prior term, see eq. (20), with η = 5 • 10 -3 δ D . Blue triangles compare the inferred fields with (x-axis) and without (y-axis) prior term. Red circles compare the inferred fields with regularization (x-axis) with the expected ones (yaxis): X DD

max 1 -

 1 αλ < , 1 -αλ > where λ >/< are the largest/smallest Eigenvalue. This equation can be solved by equating the two expression in the max:

(

  {X a } D a=1 = {{h i } N i=1 , {J ij } N i<j=1 }).As synthetic model we chose a diluted disorder Ising model on a Erdős-Rényi random graph with average connectivity c = 4:

l 4 l

 4 -l * δP P ( δP ) = (-35 ± 3.0) • 10 -4 D 2B + l -l * δP P ( δP ) = (-1.2 ± 3.0) • 10 --l DD δP P ( δP ) = (-0.3 ± 0.6) • 10 -4

ACKNOWLEDGMENTS U.F. thanks S. Deny, G. Gardella, O. Marre, R. Monasson, T. Mora, T. Obuchi and B. Telenczuk for useful discussions, S. Deny and O. Marre for the retina dataset and A. Destexhe for hosting at the European Institute for Theoretical Neuroscience. This research was supported by a grant from the Human Brain Project (HBP CLAP)

Appendix A: The existence and uniqueness of the inference solution

Because the matrix H[X], see eq. ( 9), measures the concavity of maximization problem (3), its positiveness guarantees the existence and uniqueness of the inference solution X * . As H[X] can be expressed as the covariance matrix χ[X], see eq. [START_REF] Tkacik | Searching for collective behaviour in a network of real neurons[END_REF], it is by construction nonnegative, but it could have zero modes. A zero mode in a covariance matrix identifies a linear combination of the observables, the corresponding Eigenvector, that does not fluctuate within the model probability distribution. This could happens either for all the values of the parameter X when some of the functions Σ(σ) are linearly dependent either at the solution X * when some of the fields diverge quenching (part of) the system [START_REF] Ackley | A learning algorithm for boltzmann machines[END_REF][START_REF] Barton | Ising models for neural activity inferred via selective cluster expansion: structural and coding properties[END_REF]. This last case usually happens when the dataset suffers of unsersampling.

As an example, consider a dataset of two non constant spins σ 1 (b) and σ 2 (b) that within the dataset are never active together, so that σ 1 (b)σ 2 (b) = 0 for b = 1, . . . , B. If we take only Σ 1 (σ) = σ 1 and Σ 2 (σ) = σ 2 as observables, the inference problem is well posed. However if we include, for example, Σ As introduced in sect. III A, for X ≈ X * the deterministic dynamics of the VG algorithm is exactly solvable upon projecting the fields on the χ[X * ] Eigenvectors. Along a µ-Eigenspace the convergence of the VG algorithm is not uniform and scales with the corresponding Eigenvalue λ µ as (1 -αλ µ ) t . Indeed by tuning α to speed up some particular direction, the others can suffers As introduced in sect. III B, the stochastic dynamics of the fields X is ruled by the discrete-time master equation:

where the transition rates depend on the distribution of ∇l MC X .

As Q MC X are the average of the observables Σ, for large M we expect them to be almost Gaussian distributed with a covariance equal to the Σ covariance, namely χ[X], divided by the number of MC measurements:

Consequently, as ∇l MC X ≡ B P -Q MC X we have: 26) is valid and χ[X ≈ X * ] ≈ χ , we can replace the gradient mean by the derivative of the approximated log-likelihood,

to finally obtain:

The transition rates W X→X are the probability to measure a value of l MC X such that the next fields value in the dynamics is X . From eqs. ( 21) and ( 24) it follows:

By asking P t (X) to be invariant under the evolution [START_REF] Ricci-Tersenghi | The Bethe approximation for solving the inverse Ising problem: a comparison with other inference methods[END_REF] we can obtain the stationary distribution P ∞ (X):

where δ D is the identity matrix in dimension D. Here the typical fluctuations of X around X * must consistently verify the approximation X ≈ X * : (X -X * ) 2 P∞ should be small enough to allow the expansion [START_REF] Tanaka | Mean-field theory of Boltzmann machine learning[END_REF].

We can obtain conditions on the algorithm convergence by requiring P ∞ (X) to be a properly defined probability distribution, namely to be integrable. Remarkably by asking P ∞ (X) to have a positive covariance, we re-obtain the upper-bounds for the learning rate α: αλ µ < 2, for all µ, for the VG and α < 2 for the DD.

In the stationary regime, the fluctuating fields X will induce a second source of noise in the actual distribution of Q MC . By inserting the approximation (C2) in the distribution (C1) and then by averaging over P ∞ (X), we obtain:

P ∞ (X), see eq.s ( 29) and (30), allows us to compute the expected deviation of the log-likelihood from l[X * ]. By averaging δl[X], see eq. ( 26), we obtain:

for the DD algorithm and

for the VG (we do not report the slightly involved expression of the variance, which also scales as 1/M ). Both estimations are indeed biased to lower values with large fluctuations (of the order of the bias itself). For the DD algorithm the bias depends just on D, α and M and consequently it is data independent. For the VG, instead, it depends strongly on the spectrum of χ and for large To better understand the consequences of over-fitting we consider now the case where the system that generates the data is of MaxEnt form with some unknown true fields X. By ideally sampling the distribution infinitely many times we can access to the true means of the conjugated observable P and from these compute the true log-likelihood:

We like to compare the true log-likelihood with that obtained by the exact inference, the one leading to X * and with that obtained by the DD algorithm which samples from the posterior. By sampling B times from the true model distribution we can generate synthetic dataset and obtain empirical estimates P of P. Through the central limit theorem we can approximate the distribution of the expected deviation of the observable means:

where χ = χ[ X] is the susceptibility matrix of true model. The exact inference of the fields perfectly reproducing P overestimates the log-likelihood, in fact:

where we approximate ln Z[ X + δX] up to the second order. Through the expression (D2) we can approximate the distribution of l * δP over many realization of the synthetic experiment:

where as before, D is the dimension of the fields vector and we discard terms of order B -2 . l * is on average positively biased by a factor D 2B .

- In the case of the DD algorithm (α = 1 and M = B) we have to substitute the maximization over X with an integration over the stationary fields distribution (30), which in this case will read:

where the mean equals value of δX after the maximization in the calculation of l * δP , see eq. (D3). For the DD algorithm we obtain:

and again through (D2) we obtain:

where again we discard terms of order B -2 . Coherently, the integration over the posterior distribution, that prevents to exactly maximize the log-likelihood, cancels the bias in the average of (D7).

To test these results we perform an analysis on a pairwise Ising model of N = 10 units, thus restricting the observables list to single and pairwise correlations ({Σ a } D a=1 = {{σ i } N i=1 , {σ i σ j } N i<j=1 }) and consistently the fields to biases and pairwise interactions


Clément Dumont

Giulia Baciarello

Pierre-Olivier Bosset

Pernelle Lavaud

Emeline Colomba

Christophe Massard

Yohann Loriot

Laurence Albiges

Pierre Blanchard

Alberto Bossi

Soazig Nenan

Karim Fizazi
email: karim.fizazi@gustaveroussy.fr

Prof Karim Fizazi

Long-term Castration-related Outcomes in Patients With High-risk Localized Prostate Cancer Treated With Androgen Deprivation Therapy With or Without Docetaxel and Estramustine in the UNICANCER GETUG-12 Trial

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

Introduction

Prostate cancer is the most common cancer in developed countries, with a majority of patients presenting with localized disease at diagnosis [START_REF] Horwich | Prostate cancer: ESMO Consensus Conference Guidelines[END_REF] . Treatment modalities for localized prostate cancer are based on risk factors predictive for cancer recurrence and death from the disease, mainly the TNM stage, the Gleason score, and serum PSA level. In patients with high-risk localized prostate cancer (HRLPC), the current standard of care is mostly based on prostate-focused external beam radiation therapy (EBRT) combined with androgen deprivation therapy (ADT) with LH-RH agonists for 18 to 36 months [START_REF] Nabid | Duration of Androgen Deprivation Therapy in High-risk Prostate Cancer: A Randomized Phase III Trial[END_REF][START_REF] Tosco | Systematic Review of Systemic Therapies and Therapeutic Combinations with Local Treatments for High-risk Localized Prostate Cancer[END_REF][START_REF] Gillessen | Management of Patients with Advanced Prostate Cancer: The Report of the Advanced Prostate Cancer Consensus Conference APCCC 2017[END_REF] .

Chemotherapy with docetaxel given concomitantly with ADT is a time-honored standard for patients with metastatic castration-resistant prostate cancer [START_REF] Tannock | Docetaxel plus prednisone or mitoxantrone plus prednisone for advanced prostate cancer[END_REF] . More recently, this treatment has also been shown to improve overall survival when administered in the castration-naive setting [START_REF] Sweeney | Chemohormonal Therapy in Metastatic Hormone-Sensitive Prostate Cancer[END_REF][START_REF] James | Addition of docetaxel, zoledronic acid, or both to first-line long-term hormone therapy in prostate cancer (STAMPEDE): survival results from an adaptive, multiarm, multistage, platform randomised controlled trial[END_REF] . Whether combining chemotherapy with ADT even earlier in the natural history of prostate cancer might benefit the patients, it has been a subject of inquiry for almost two decades [START_REF] Vale | Addition of docetaxel or bisphosphonates to standard of care in men with localised or metastatic, hormone-sensitive prostate cancer: a systematic review and metaanalyses of aggregate data[END_REF] .

The UNICANCER GETUG 12 trial was the first phase III trial to assess the role of docetaxel-based chemotherapy in men with localized prostate cancer: 413 patients with previously untreated HRLPC were randomized to receive ADT with or without concomitant chemotherapy with docetaxel plus estramustine (DE) before local treatment, which turned out to be EBRT in 87% of patients [START_REF] Fizazi | A phase III trial of docetaxel-estramustine in high-risk localised prostate cancer: A planned analysis of response, toxicity and quality of life in the GETUG 12 trial[END_REF] . The addition of chemotherapy resulted in significantly improved relapse-free survival [START_REF] Fizazi | Androgen deprivation therapy plus docetaxel and estramustine versus androgen deprivation therapy alone for high-risk localised prostate cancer (GETUG 12): a phase 3 randomised controlled trial[END_REF] , the primary endpoint of the trial. With a median follow-up of 12 years, this trial showed that DE reduces biochemical failures, but also improves clinical or radiological relapse-free survival with events including metastases, local relapses, and deaths [START_REF] Fizazi | Updated results of GETUG-12, a phase III trial of docetaxel-based chemotherapy in high-risk localized prostate cancer, with a 12-year follow-up[END_REF] . Since GETUG-12 primary analysis read out in 2015, other trials randomly testing docetaxel in men with HRLPC have also reported early relapse-free survival, with some showing relapse-free survival improvement [START_REF] Tosco | Systematic Review of Systemic Therapies and Therapeutic Combinations with Local Treatments for High-risk Localized Prostate Cancer[END_REF][START_REF] James | Addition of docetaxel, zoledronic acid, or both to first-line long-term hormone therapy in prostate cancer (STAMPEDE): survival results from an adaptive, multiarm, multistage, platform randomised controlled trial[END_REF][START_REF] Schweizer | Adjuvant Leuprolide With or Without Docetaxel in Patients With High-Risk Prostate Cancer After Radical Prostatectomy (TAX-3501)[END_REF][START_REF] Lin | VA CSP#553 Chemotherapy after prostatectomy (CAP) for hisgh-risk prostate carcinoma: a phase III randomized study[END_REF][START_REF] Ahlgren | Docetaxel Versus Surveillance After Radical Prostatectomy for High-risk Prostate Cancer: Results from the Prospective Randomised, Open-label Phase 3 Scandinavian Prostate Cancer Group 12 Trial[END_REF][START_REF] Rosenthal | Effect of Chemotherapy With Docetaxel With Androgen Suppression and Radiotherapy for Localized High-Risk Prostate Cancer: The Randomized Phase III NRG Oncology RTOG 0521 Trial[END_REF] . The benefit of adjuvant/neoadjuvant docetaxel on relapse-free survival was further confirmed by a meta-analysis including the first four trials with available data [START_REF] Vale | Addition of docetaxel or bisphosphonates to standard of care in men with localised or metastatic, hormone-sensitive prostate cancer: a systematic review and metaanalyses of aggregate data[END_REF] . In the GETUG 12 trial, the on-treatment toxicity of chemotherapy was manageable, with only a 2% febrile neutropenia rate, and its negative impact on quality of life was only transient, with patient-reported outcomes being similar in both arms one year after the end of DE treatment [START_REF] Fizazi | A phase III trial of docetaxel-estramustine in high-risk localised prostate cancer: A planned analysis of response, toxicity and quality of life in the GETUG 12 trial[END_REF] . These observations are supported by similar data obtained in the metastatic setting [START_REF] Morgans | Quality of Life During Treatment With Chemohormonal Therapy: Analysis of E3805 Chemohormonal Androgen Ablation Randomized Trial in Prostate Cancer[END_REF] .

Little is known about the very long-term impact of early taxane-based chemotherapy when combined with ADT. Here, we present the results of a longitudinal study of patients from Institute Gustave Roussy treated in the UNICANCER GETUG 12, focusing on recovery from castration and its associated functional and cardiovascular outcomes.

Methods

Patients and treatment

The detailed inclusion criteria and treatment plan for the GETUG 12 trial have been previously described [START_REF] Fizazi | A phase III trial of docetaxel-estramustine in high-risk localised prostate cancer: A planned analysis of response, toxicity and quality of life in the GETUG 12 trial[END_REF] . Briefly, patients had treatment-naive, histologically proven prostate cancer, no evidence of distant metastases on both abdomino-pelvic CT and bone scan performed during the last 6 months, and one or more adverse risk features among: clinical stage T3 or T4, serum PSA level >20 ng/mL, Gleason score ≥8, histologically proven regional lymph node involvement. All patients underwent surgical lymph node staging during the 12 weeks before randomization. All patients were treated with subcutaneous goserelin 10.8 mg every 3 months for 3 years, concomitantly with an androgen receptor inhibitor for the first 3 weeks. For patients allocated to the experimental arm (ADT+DE arm), chemotherapy with intravenous docetaxel 70 mg/m² on day 2 associated with oral prednisone on day 1 to 3 and oral estramustine 10 mg/kg/day on day 1 to 5 were administered every 3 weeks for 4 cycles before local treatment, with concomitant thromboprophylaxis. In both arms, local treatment was performed 3 months after the start of ADT and could consist of either EBRT or prostatectomy for patients with node-negative disease, and either EBRT or no local treatment for patients with node-positive disease. Radiotherapy procedures have been previously described [START_REF] Fizazi | A phase III trial of docetaxel-estramustine in high-risk localised prostate cancer: A planned analysis of response, toxicity and quality of life in the GETUG 12 trial[END_REF] .

Data collection

For the present study, the subgroup of patients enrolled in the GETUG 12 trial at the Institut Gustave Roussy (Villejuif, France) were monitored twice yearly as per protocol after ADT discontinuation. Supplementary clinical data regarding body weight, libido, and quality of erection were collected prospectively during that follow-up. Libido was characterized as either absent, decreased, or normal.

Quality of erection was characterized as no erection, erection present but not permitting intercourse (patient complaining or not complaining), erection permitting intercourse with medical assistance, or erection permitting intercourse without medical assistance. Serum testosterone levels (STLs) were collected prospectively. Data was collected annually for up to 10 years after randomization (i.e. 7 years post-ADT). Patients were excluded from the analyses in case of death, prostate cancer recurrence leading to ADT resumption, or any other illness that may interfere with the interpretation of the aforementioned parameters. Baseline data was collected retrospectively using the local computerized individual patient files. Cardiovascular events notified in these files were collected.

Outcomes description and statistical analysis

The distribution of individual changes in body weight between randomization and ADT discontinuation in both treatment arms was compared using a Students's t-test. The percentage of patients with a 10% or more weight gain at the end of ADT was compared using a chi2 test. The timely evolution of the mean change in body weight in each treatment arm was described.

Discreet STL values were dichotomized as either castrate (<50 ng/dL) or supracastrate (≥50 ng/dL), and STLs ≥200 ng/dL were considered as normal [START_REF] Takei | Duration of androgen deprivation therapy and nadir of testosterone at 20 ng/dL predict testosterone recovery to supracastrate level in prostate cancer patients who received external beam radiotherapy[END_REF] . The proportion of patients with supracastrate and normal STLs was described at the end of ADT, 6 months post-ADT, and then yearly. The proportion evaluable patients recovering to normal STLs at any time after ADT discontinuation in both treatment arms was compared using a chi2 test.

Libido was analyzed as a binary parameter, either normal or absent/decreased. The proportion of evaluable patients reporting normal libido at least once at any time after ADT discontinuation in both treatment arms was compared using a chi2 test.

Quality of erection was analyzed as a binary parameter, either erection permitting intercourse without medical assistance or other. The proportion of evaluable patients reporting erections allowing for intercourse without medical help any time after ADT discontinuation in both treatment arms was compared using a chi2 test. For homogeneity concerns, patients treated with radical prostatectomy were excluded from this analysis.

Results

Patients and treatment

Between November 2002 and November 2006, 81 patients were included in the GETUG 12 trial at the Institute Gustave Roussy; 2 patients had died and one patient was lost to follow-up before completion of the 3 years of ADT. Of the remaining 78 patients, 36 patients had been randomized to the ADT arm and 42 patients to the ADT+DE arm. Baseline characteristics as well as the nature of local treatment received by the patients are displayed in Table 1.

At data cutoff on March 15, 2014, the median follow-up after ADT discontinuation was 5.9 years.

Patients were censored from the analyses before data cutoff for the following reasons: prostate cancer recurrence needing ADT resumption (n=6) or salvage radiation therapy (n=1), other medical condition (n=2), death (n=8), or lost to follow-up (n=4).

Changes in body weight

At the end of ADT, 21% and 15% of evaluable patients presented a >10% weight gain in the ADT arm and the ADT+DE arm, respectively (p=0.51). There was no difference in the distribution of changes in body weight from baseline in both treatment arms (t-test: p=0.47; Figure 1A). The mean weight variation at the end of ADT was similar in both treatment arms (+3.9% and +3.5%, respectively), with a similar trend toward timely normalization in patients with prolonged follow-up (Figure 1B).

Testosterone recovery

As expected, STLs were constantly at castrate levels at the end of ADT in all evaluable patients. The proportion of patients with supracastrate STLs increased to 41%, 81%, and 91% of evaluable patients after 6 months, 1 year, and 3 years of ADT discontinuation, respectively (Figure 2A), and the proportion of patients with normal STLs were 29%, 73%, and 81%, respectively (Figure 2B). In 72 evaluable patients with at least one available STL value after ADT discontinuation, patients recovering to normal STLs at any time point after ADT discontinuation was 80% in the ADT arm and 78% in the ADT+DE arm (p=0.87).

Functional outcomes

Similarly, the proportion of patients reporting normal libido and erections allowing for intercourse without medical help increased during the first 4 years following ADT discontinuation (Figure 3). Among evaluable patients, 44% of patients in the ADT arm and 59% in the ADT+DE arm reported normal libido at some point after ADT discontinuation (p=0.20), and 37% and 47% reported erections allowing for intercourse without medical help, respectively (p=0.37).

Cardiovascular morbidity

As of data cutoff, an arterial thrombotic event had occurred in 6 patients and a venous thrombotic event had occurred in 4 patients (summarized in Table 2). Although this number of events was too small to allow for statistical analysis, no gross excess of cardiovascular events was apparent in one or the other treatment arm.

Discussion

This study provides the first insight on long-term castration-related outcomes in patients treated with upfront docetaxel-based chemotherapy and long-course ADT for HRLPC. The key observation derived from our data is that neither changes in body weight, testosterone recovery, libido or quality of erection seemed to be negatively impacted by chemotherapy. The long follow-up and the randomized nature of this study strengthen these findings.

In our study, restoration of STLs to normal levels was achieved in 78% of evaluable patients within 2 years of ADT discontinuation, then plateaued. Despite methodological discrepancies, the observations seem to be coherent with the rare previous reports of testosterone recovery after longcourse adjuvant/neoadjuvant ADT [START_REF] Yoon | Testosterone Recovery After Prolonged Androgen Suppression in Patients With Prostate Cancer[END_REF][START_REF] Planas | Hormonal response recovery after long-term androgen deprivation therapy in patients with prostate cancer[END_REF] . Although the number of patients was too small to draw any definitive conclusions, baseline body mass index (BMI) seemed to negatively impact testosterone recovery, with evaluable patients with BMI <25 kg/m², 25-29.9 kg/m², and ≥30 kg/m² recovering normal STLs in 87%, 79%, and 67% of cases, respectively. Weight gain on ADT was not associated with a decreased probability of normal testosterone recovery. One limit of our study was the absence of data regarding baseline STLs, which precludes analyzing non-recovery at the individual level; indeed, previous studies have reported that pre-ADT testosterone levels are predictive of testosterone recovery [START_REF] Padula | Normalization of serum testosterone levels in patients treated with neoadjuvant hormonal therapy and three-dimensional conformal radiotherapy for prostate cancer[END_REF] . Importantly, we observed no apparent delay in testosterone recovery associated with chemotherapy. This is in line with findings from the TAX-3501 trial, in which addition of docetaxel to leuprolide after radical prostatectomy (either as adjuvant therapy or differed until progression) did not impact testosterone recovery [START_REF] Schweizer | Adjuvant Leuprolide With or Without Docetaxel in Patients With High-Risk Prostate Cancer After Radical Prostatectomy (TAX-3501)[END_REF] . An important difference between the two studies is that, in the TAX-3501 trial, ADT was pursued for only 18 months when several reports have concluded that a duration of treatment longer than 18 or 24 months negatively impacts testosterone recovery in patients treated with adjuvant ADT [START_REF] Takei | Duration of androgen deprivation therapy and nadir of testosterone at 20 ng/dL predict testosterone recovery to supracastrate level in prostate cancer patients who received external beam radiotherapy[END_REF][START_REF] Yoon | Testosterone Recovery After Prolonged Androgen Suppression in Patients With Prostate Cancer[END_REF][START_REF] Nam | Factors associated with testosterone recovery after androgen deprivation therapy in patients with prostate cancer[END_REF] . Since 18 months to 3 years of adjuvant ADT are now the standard of care in HRLPC patients treated with EBRT [START_REF] Tosco | Systematic Review of Systemic Therapies and Therapeutic Combinations with Local Treatments for High-risk Localized Prostate Cancer[END_REF][START_REF] Gillessen | Management of Patients with Advanced Prostate Cancer: The Report of the Advanced Prostate Cancer Consensus Conference APCCC 2017[END_REF] , we felt that confirmation of those findings from TAX-3501 was needed. Other differences include the nature of local treatment, which was radical prostatectomy, and length of follow-up, which was only 3.4 years from treatment initiation in the report from TAX-3501. Finally, a major specificity of the GETUG 12 trial was the addition of estramustine to docetaxel. Estramustine has long been described as having endocrine effects including gonadotropic suppression [START_REF] Bishop | Plasma hormone levels in patients with prostatic carcinoma treated with diethylstilboestrol and estramustine[END_REF][START_REF] Morales | Clinical relevance of plasma testosterone and prolactin changes in advanced cancer of prostate treated with diethylstilbestrol or estramustine phosphate[END_REF] , and the long-term endocrine toxicity of combined LH-RH agonists and estramustine had been a concern. In our subset of patients from the GETUG-12 trial we observed no excess in long-term endocrine toxicity in estramustine-treated patients. The similarity in the rates and timeliness in testosterone recovery in both arms also does not support the hypothesis that the benefit associated with neoadjuvant docetaxel and estramustine in the GETUG 12 trial was due to deeper or prolonged castration.

Our study is the first to report on long-term functional outcomes in patients treated with neoadjuvant chemotherapy in combination with ADT for HRLPC. The drastic on-treatment impact of ADT on sexuality-related patient-reported outcomes (PROs) is well-described [START_REF] Potosky | Quality-of-life outcomes after primary androgen deprivation therapy: results from the Prostate Cancer Outcomes Study[END_REF][START_REF] Gay | Neoadjuvant Androgen Deprivation Therapy Leads to Immediate Impairment of Vitality/Hormonal and Sexual Quality of Life: Results of a Multicenter Prospective Study[END_REF] . In accordance with previous reports [START_REF] Marumo | Erectile function and nocturnal penile tumescence in patients with prostate cancer undergoing luteinizing hormone-releasing hormone agonist therapy[END_REF] , the castrate STLs observed in all patients at the end of ADT in our study were associated with dramatic decreased libido, and virtually no patient reported the possibility of sexual intercourse without medical assistance. Few studies have focused on such functional outcomes beyond treatment discontinuation. For patients undergoing short-term ADT and radiation therapy, conflicting reports exist as to the reversibility of ADT-induced changes in sexual activity and sexual functioning within 2 years of treatment discontinuation [START_REF] Bolla | Short Androgen Suppression and Radiation Dose Escalation for Intermediateand High-Risk Localized Prostate Cancer: Results of EORTC Trial 22991[END_REF][START_REF] Gay | External Beam Radiation Therapy or Brachytherapy With or Without Shortcourse Neoadjuvant Androgen Deprivation Therapy: Results of a Multicenter, Prospective Study of Quality of Life[END_REF] ; this may be the reflection of methodological discrepancies in PRO measurement. In any case, we feel that a 2-year follow-up may be too short to evaluate such outcomes since a decline in erectile function may be observed even later in these patients [START_REF] Daly | The effect of short term neo-adjuvant androgen deprivation on erectile function in patients treated with external beam radiotherapy for localised prostate cancer: An analysis of the 4-versus 8-month randomised trial (Irish Clinical Oncology Research Group 97-01)[END_REF] . Data regarding recovery after long-term ADT is even more scarce [START_REF] Wilke | Testosterone and erectile function recovery after radiotherapy and longterm androgen deprivation with luteinizing hormone-releasing hormone agonists[END_REF] . We did not use validated questionnaires regarding libido and quality of erection, which is a limitation of our study. Nevertheless, we were able to observe that both symptoms significantly improved over time, with about 40% of patients regaining the possibility of intercourse without medical assistance within 2 to 3 years of ADT discontinuation. We did not find any association between baseline erectile dysfunction and post-ADT functional outcomes, possibly because of a lack of statistical power and of missing data. The frequent discrepancy between normalized STLs and abnormal functional outcomes in most patients pleads for non-endocrine causes. Indeed, most men in GETUG-12 also received radiotherapy to the prostate, this may have impaired their erectile function.

Cardiovascular morbidity and mortality are chief concerns in prostate cancer patients treated with long-term ADT [START_REF] D'amico | Influence of androgen suppression therapy for prostate cancer on the frequency and timing of fatal myocardial infarctions[END_REF][START_REF] Keating | Diabetes and cardiovascular disease during androgen deprivation therapy: observational study of veterans with prostate cancer[END_REF][START_REF] Zapatero | Late Radiation and Cardiovascular Adverse Effects After Androgen Deprivation and High-Dose Radiation Therapy in Prostate Cancer: Results From the DART 01/05 Randomized Phase 3 Trial[END_REF] . Nevertheless, the incidence of cardiovascular or thrombotic events was low in GETUG-12 and, in the subgroup of patients analyzed in the present paper, we did not observe any gross increase in those randomized in the chemotherapy arm. This is reassuring regarding the longterm toxicity of neoadjuvant docetaxel, especially when considering that in the GETUG 12 trial chemotherapy also included estramustine, which is a known pro-thrombotic drug [START_REF] Lubiniecki | Thromboembolic events with estramustine phosphate-based chemotherapy in patients with hormone-refractory prostate carcinoma: results of a meta-analysis[END_REF] . Also, no excess in second cancers was observed in the recently reported long-term analysis of the trial [START_REF] Fizazi | Updated results of GETUG-12, a phase III trial of docetaxel-based chemotherapy in high-risk localized prostate cancer, with a 12-year follow-up[END_REF] .

The benefit of docetaxel-based chemotherapy in HRLPC is suggested by a building corpus of evidence, although not yet mature enough to support its implementation in routine practice in most countries [START_REF] Vale | Addition of docetaxel or bisphosphonates to standard of care in men with localised or metastatic, hormone-sensitive prostate cancer: a systematic review and metaanalyses of aggregate data[END_REF][START_REF] Tosco | Systematic Review of Systemic Therapies and Therapeutic Combinations with Local Treatments for High-risk Localized Prostate Cancer[END_REF][START_REF] Rosenthal | Effect of Chemotherapy With Docetaxel With Androgen Suppression and Radiotherapy for Localized High-Risk Prostate Cancer: The Randomized Phase III NRG Oncology RTOG 0521 Trial[END_REF] . While the long-term results of pioneer studies such as the GETUG-12 trial are made available [START_REF] Fizazi | Updated results of GETUG-12, a phase III trial of docetaxel-based chemotherapy in high-risk localized prostate cancer, with a 12-year follow-up[END_REF][START_REF] Rosenthal | Effect of Chemotherapy With Docetaxel With Androgen Suppression and Radiotherapy for Localized High-Risk Prostate Cancer: The Randomized Phase III NRG Oncology RTOG 0521 Trial[END_REF] , other studies are ongoing such as the PEACE-2/GETUG-AFU 23 trial (NCT01952223) which explores the benefit of neoadjuvant cabazitaxel, a taxane active in metastatic patients pretreated with docetaxel [START_REF] De Bono | Prednisone plus cabazitaxel or mitoxantrone for metastatic castrationresistant prostate cancer progressing after docetaxel treatment: a randomised open-label trial[END_REF] , in HRLPC patients treated with EBRT with or without pelvic irradiation.

The data we derived from a subset of the GETUG 12 trial population is reassuring as of the long-term endocrine or cardiovascular side effects associated with additional chemotherapy in this setting and encourages pursuing such trials. Nevertheless, because the present study is a subgroup analysis, results are only informative and need to be interpret cautiously. In our opinion, any future studies of early chemotherapy in HRLPC patients large enough to address these issues, such as the PEACE-2 trial, should include a prospective survey of the long-term biochemical and functional recovery from castration, as well as of metabolic and cardiovascular morbidity.

Conclusion

Additional neoadjuvant chemotherapy with docetaxel and estramustine did not negatively impact neither testosterone recovery nor long-term functional outcomes, changes in body weight or the incidence of cardiovascular events in the first years following treatment discontinuation after longcourse androgen deprivation therapy for high-risk localized prostate cancer.

Figure 1 -Changes in body weight

Table 1 -Baseline characteristics for patients alive at the end of ADT

 1 A. Distribution of changes in body weight at the time of ADT discontinuation (ADT arm: 29 evaluable patients; ADT+DE arm: 41 evaluable patients) B. Timely evolution of mean weight variation from baseline after ADT discontinuation (numbers of evaluable patients in each treatment arm at each time point are shown) DE: docetaxel and estramustine ; BMI: body mass index ; PSA: prostate-specific antigen ; EBRT: external beam radiotherapy * Atheromatosis include coronary artery disease, carotid artery disease, and lower limbs artery disease

	50%
	40%	ADT
	30%	ADT+CT
	20%
	10%
	0%
		≥10% loss	5-10% loss	Stable	5-10% gain	≥10% gain

Acknowledgements

The authors would like to thank Dr Lilian Amrein for her writing assistance and liguistical help.


Olivier Plateker
email: olivier.plateker@gmail.com

Thomas Daix

MD Cyrille Boulogne

Virginie Rousselle

MD Philippe Vignon

PhD Jean Porterie

Post-kyphoplasty cement embolism migrating to the peritoneum through the right ventricle

Keywords: kyphoplasty, embolism, pericardium

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

2

A 58-year-old woman was admitted to the Emergency Department for dyspnea, nonproductive cough and acute left shoulder pain, radiating to the left lower anterior chest, and occurring two days after a kyphoplasty of the first and third lumbar vertebrae for osteoporotic fractures. Physical examination revealed a non-productive cough but no sign of respiratory distress or cyanosis, normal cardiopulmonary auscultation and normal heart rate and blood pressure. ECG disclosed a sinus rhythm with right bundle branch block and related nonspecific repolarization changes. Contrast-enhanced CT-scan excluded pulmonary embolism but disclosed a high-attenuation 6-cm foreign body traversing the right ventricle and extending into the pericardium, without extravasation of contrast but with mild pericardial effusion (Fig A andB). Transthoracic echocardiography confirmed the presence of the foreign body, up to 7 mm thick, within the pericardial space in conjunction with mild pericardial effusion, but neither signs of tamponade nor regional wall motion abnormality ( Kyphoplasty is a per-cutaneous technique mainly indicated for the treatment of osteoporotic or metastatic vertebral fracture. Although minimally invasive, the procedure may lead to cement contamination of the venous system and spinal canal. In a recent meta-analysis incorporating 22 studies, Zhan et al. [START_REF] Zhan | Risk Factors for Cement Leakage After Vertebroplasty or Kyphoplasty: A Meta-Analysis of Published Evidence[END_REF] reported an incidence of cement leakage of 60% after percutaneous vertebroplasty and 18% after percutaneous balloon kyphoplasty. Intracardiac cement embolism has been scarcely reported in the literature as opposed to pulmonary embolism. In a retrospective study including 1548 patients 2 intracardiac cement embolism was diagnosed in 72 patients (3.9%), but was symptomatic in only 6 of them and always associated with pulmonary embolism. As far as we know, intraperitoneal migration of cement embolism after cardiac perforation without tamponade or pulmonary embolism has not yet been described.

We report the first case of intraperitoneal migration of cement embolism through right ventricular and diaphragm perforation after kyphoplasty, without associated tamponade or pulmonary embolism.

  Fig C). Review of the fluoroscopic images acquired at the end of the kyphoplasty confirmed the presence of a cement embolism within a vertebral vein arising at L3 (Fig D). Twenty-four hours later, a new CT-scan depicted migration of the foreign body toward the peritoneum through the diaphragm (Fig D), and a truncated aspect of the cement link from the lumbar vertebra (Fig E). The cement embolism was extracted from the peritoneal cavity through a sub-xiphoid approach (Fig F). Post-operative course was uneventful and the patient was discharged at home on day 10.