Said Obakrim
email: said.obakrim@univ-rennes1.fr

Pierre Ailliot
email: pierre.ailliot@univ-brest.fr

Valérie Monbet
email: valerie.monbet@univ-rennes1.fr

Nicolas Raillard
email: nicolas.raillard@ifremer.fr

EM ALGORITHM FOR GENERALIZED RIDGE REGRESSION WITH SPATIAL COVARIATES

Keywords: Generalized Ridge, EM algorithm, Spatial covariates, Matérn, Conditional Autoregressive

The generalized Ridge penalty is a powerful tool for dealing with overfitting and for high-dimensional regressions. The generalized Ridge regression can be derived as the mean of a posterior distribution with a Normal prior and a given covariance matrix. The covariance matrix controls the structure of the coefficients, which depends on the particular application. For example, it is appropriate to assume that the coefficients have a spatial structure in spatial applications. This study proposes an expectation-maximization algorithm for estimating generalized Ridge parameters whose covariance structure depends on specific parameters. We focus on three cases: diagonal (when the covariance matrix is diagonal with constant elements), Matérn, and conditional autoregressive covariances. A simulation study is conducted to evaluate the performance of the proposed method, and then the method is applied to predict ocean wave heights using wind conditions.

Introduction

Consider an experiment where we have the data {y, X}, of n observations of a continuous variable Y and n × d matrix of covariates X. Suppose that Y is related to X via a linear model

Y = Xβ + , (1)
where β are model coefficients and ∼ N (0, σ 2 ) is the model error. We suppose that the intercept is either included in β (so that the first column of X is a vector of 1) or that Y and X are centered. The least squares estimates are the best linear unbiased estimates of the parameters β. However, in the case of multicollinearity or high-dimensionality, penalized linear regression methods, like Ridge regression, are needed to control the variance. Ridge estimator of the problem (1) is βRidge

λ = arg min β -(β, σ 2 ) + λ β 2 (2)
where λ is the regularization parameter and (β, σ 2 ) is the log-likelihood of the model (1). High values of λ permit to reduce the variance and increase the bias of the model. A good model should have a trade-off between variance and bias [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF]. In order to find a trade-off between bias and variance, the hyperparameter λ needs to be selected. [START_REF] Boonstra | A small-sample choice of the tuning parameter in ridge regression[END_REF] classified methods for selecting λ into goodness-of-fit-based and likelihoodbased methods. Goodness-of-fit-based methods define a goodness of fit criterion (such as the mean squared error) and minimize it in terms of λ. The most common goodness-of-fit-based method is the k-fold cross-validation which consists of partitioning observations into k groups and estimating β k times for each λ leaving out one group. For each λ, a goodness of fit score is calculated, and λ with the maximum score value is chosen. The typical choice of k is 5 and 10, while setting k = n leads to leave-one-out cross-validation (LOOCV). LOOCV leads to a better estimation of λ; however, it is computationally expensive given that it requires fitting the model n times [START_REF] Patil | Uniform consistency of cross-validation estimators for highdimensional ridge regression[END_REF]. Generalized cross-validation (GCV) [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF]) is an approximation of LOOCV that does not require fitting n models. GCV uses a weighted version of the predicted residual error sum of squares (PRESS) statistic [START_REF] Allen | The relationship between variable selection and data agumentation and a method for prediction[END_REF] as a goodness of fit criterion. One of the problems with goodness-of-fit-based methods is the selection of the grid of λ, which influences the estimation.

Assuming that Y |β ∼ N (Xβ, σ 2 I n ), Ridge regression can be derived as the mean of a posterior distribution with the prior β ∼ N (0 d , σ 2 λ -1 I d ) [START_REF] Van Wieringen | Lecture notes on ridge regression[END_REF] and as in Bayesian hierarchical linear regression, likelihood-based methods maximize the likelihood with respect to σ 2 and λ using for instance an iterative method [START_REF] Lee | Hierarchical generalized linear models[END_REF]. Unlike goodness-of-fit-based methods, the advantage of likelihood-based approaches is, on the one hand, that they do not require grid selection for the regularization parameters. On the other hand, likelihood-based methods can be generalized to consider any form of prior for the coefficients β. In some applications, the regression coefficients can be penalized differently, or a joint penalization of the coefficients is required. For example, in spatial statistics, where predictors have a spatial structure, it is reasonable to suppose that coefficients have a spatial structure. To do that, the generalized Ridge (van Wieringen, 2015) can be used. Generalized Ridge extends the equation ( 2) by replacing the term λ β 2 to β T ∆β, where ∆ is called the penalty matrix. In general, ∆ depends on some regularization parameters (see, e.g., [START_REF] Goeman | Autocorrelated logistic ridge regression for prediction based on proteomics spectra[END_REF] and [START_REF] Hemmerle | An explicit solution for generalized ridge regression[END_REF]); however, when the number of the regularization parameters is greater than 1, goodness-of-fit-based methods struggle with the problem of combinatorial explosion. Generalized Ridge in the hierarchical linear model framework, is equivalent to suppose that β ∼ N (0 d , Σ θ ) where Σ θ is a covariance matrix that depends on some parameters θ. Note that Σ θ corresponds to the inverse of the penalty matrix ∆. The classical Ridge is a special case of this model when the covariance matrix Σ θ is diagonal, and θ is the usual regularization parameter λ.

Considering β as a hidden variable, [START_REF] Bishop | Pattern recognition and machine learning[END_REF] proposed an expectation-maximization (EM) algorithm to find the maximum likelihood estimation (MLE) of parameters of a Bayesian linear regression model. The EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF] is a method for estimating the parameters of a model with hidden variables. The EM algorithm alternates between two steps: the expectation and maximization steps. The E-step calculates the conditional expectation of the log-likelihood given the observations and current parameters. In the M-step, the parameters are estimated by maximizing the conditional expectation of the log-likelihood calculated in the E-step.

In this study, we extend the algorithm in [START_REF] Bishop | Pattern recognition and machine learning[END_REF] and propose an EM algorithm to estimate the parameters of hierarchical linear regression when β ∼ N (0, Σ θ ). At first, we study the case where Σ θ is diagonal with constant elements, which corresponds to the classical Ridge in equation ( 2) and the problem studied by [START_REF] Bishop | Pattern recognition and machine learning[END_REF]. Then, we consider the case where the coefficients β have a spatial structure, especially when Σ θ is the Matérn or the conditional autoregressive (CAR) covariance. A simulation study is done to assess the performance of the method. Then, the proposed method is applied to oceanographic data where the response variable represents a wave parameter in a location in the Bay of Biscay, and X represents wind conditions over the North Atlantic [START_REF] Obakrim | Statistical modeling of the space-time relation between wind and significant wave height[END_REF].

This paper is organized as follows. The proposed method and its special cases are presented in Section 2. Then, a simulation study is conducted in Section 3 to assess the performance of the proposed method. In section 4, we apply the methodology to oceanography data. Finally, this study is concluded in Section 5.

Proposed method

As stated in the introduction, Ridge regression can be viewed as a hierarchical linear model where

β ∼ N (0 d , σ 2 λ -1 I d ).
When there is a structure on the coefficients, it is unreasonable to consider all possible covariance functions as possible candidates for β. Therefore, we suppose that the covariance of β depends on some parameters θ, so that β ∼ N (0 d , Σ θ ). This motivates using the EM algorithm to find the maximum likelihood estimation of the parameters, where the model parameters are then Θ = (σ 2 , θ). The proposed method is described in this section, and three special cases of the covariance Σ θ (the diagonal, Matérn, and CAR) are studied.

EM algorithm for generalized Ridge

Consider the linear model (1) and assume that β is a latent variable that follows a normal distribution. We define the regression model hierarchically as

β ∼ N (0 d , Σ θ ) Y | β, Θ ∼ N (Xβ, σ 2 I n ) (3)
where Θ = (σ 2 , θ). Note that for simplicity, we assume that the mean of β is zero. The EM algorithm for the case where β has a non-zero mean will be presented in the Appendix.

Given a sample y = (y 1 , ..., y n ), the complete log-likelihood is expressed as ln p(y, β; Θ) = ln p(y | β; σ 2 ) + ln p(β; θ)

= - 1 2 d ln(2π) + ln(|Σ θ |) + β T Σ -1 θ β + n ln(2π) + n ln(σ 2 ) + 1 σ 2 y -Xβ 2 (4)
Maximum likelihood estimation consists of maximizing (4) with respect to the parameters Θ. This is usually done with the Expectation-Maximization algorithm in the latent variable context. The EM algorithm alternates between the E-step and M-step. In the E-step, the expectation Q(Θ|Θ (t) ) of the complete likelihood with respect to the posterior distribution of the latent variable β and the parameters Θ (t) from the previous iteration t is calculated. In the M-step, the quantity Q(Θ|Θ (t) ) is maximized with respect to the parameters Θ.

The E-step and M-step are defined as follows

• E-step: t) ).

Q(Θ|Θ (t) ) = E(ln p(y, β; Θ) | y, Θ (
(5)

The posterior distribution of the latent variable β is a normal distribution with mean µ β|y and covariance matrix Σ β|y such that

Σ β|y = (Σ -1 θ + 1 σ 2 X T X) -1 µ β|y = (X T X + σ 2 Σ -1 θ ) -1 X T y. (6)
Note that µ β|y defined in ( 6) is a generalized Ridge estimator (see e.g. van Wieringen ( 2015)) solution of the optimization problem

µ β|y = arg min β y -Xβ 2 σ 2 + β T Σ -1 θ β (7)
Therefore,

Q(Θ|Θ (t) ) = - 1 2 ln(|Σ θ |) + Tr(Σ -1 θ E(ββ T | y, Θ (t) )) + ln(σ 2 ) + 1 σ 2 E( y -Xβ 2 | y, Θ (t) ) + C (8)
where C is a constant and

E(ββ T |y; Θ (t) ) = Σ β|y + µ β|y µ T β|y E( y -Xβ 2 |y; Θ (t) ) = y 2 -2y T Xµ β|y + Tr(X T XE(ββ T |y; Θ (t) )) (9)
• M-step:

The maximization step computes

Θ (t+1) = arg max Θ Q(Θ|Θ (t) ) (10)
which leads to the following updates of the parameters σ 2 and θ

σ 2,(t+1) = 1 n ( y 2 -2y T Xµ β|y + Tr(X T XE(ββ T |y; Θ (t) ))) θ (t+1) = arg max θ ln(|Σ -1 θ |) -Tr(Σ -1 θ E(ββ T | y, Θ (t) )) (11)

Special cases

The M-step in equation ( 11) requires the maximization of Q(Θ|Θ (t) ) over the parameters of the covariance Σ θ . In this study, we will explore three cases. First, we consider the case where Σ θ is diagonal. Then, the case where β has a spatial structure, especially when the parametric covariance is the Matérn covariance function. Finally, we consider the conditional autoregressive model (CAR).

Diagonal case

In the classical Ridge, the covariance matrix of the coefficients β is supposed to be diagonal such that

Σ θ = σ 2 β I d . (12)
The M-step of the covariance in (11) becomes

σ 2,(t+1) β = arg max σ 2 β -d ln(σ 2 β ) - 1 σ 2 β Tr(E(ββ T | y, Θ (t) )). ( 13
)
Setting the derivatives with respect to σ 2 β to zero, we obtain the M-step

σ 2,(t+1) β = Tr(E(ββ T | y, Θ (t) )) d . ( 14
)
Note that 1 σ 2 β corresponds to the regularization parameter λ in equation ( 1). As stated in the introduction, Ridge regression requires the selection of the regularization parameter. Therefore, the EM algorithm can be an alternative to cross-validation for estimating Ridge coefficients along with the regularization parameter. A comparison of the two methods (cross-validation and EM algorithm) is given in the Appendix.

Spatial covariance functions

In spatial statistics applications, one may assume that β has a spatial structure. One way to do that is to assume that β has a parametric covariance function. There are many choices of covariance functions that are widely used for Gaussian processes and kriging [START_REF] Schulz | A tutorial on gaussian process regression: Modelling, exploring, and exploiting functions[END_REF]. In this study, we focus on the stationary Matérn covariance, which has the form

K(h; φ, κ) = σ 2 β 2 κ-1 Γ(κ) h φ κ K κ h φ ( 15
)
where h is the distance between two points, Γ is the Gamma function, and K κ is the modified Bessel function [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]. The Matérn function is parameterized by the variance parameter σ 2 β , the range parameter φ, and the smoothness parameter κ. The range parameter φ controls the decay rate with distance, with larger values of φ corresponding to more strongly correlated variables, and the smoothness parameter κ controls the mean-square differentiability of the spatial process.

The M-step of the covariance of β in (12) becomes

(σ 2,(t+1) β , θ (t+1) ) = arg max σ 2 β ,θ ln(|R -1 θ |) -d ln(σ 2 β ) - 1 σ 2 β Tr(R -1 θ E(ββ T | y, Φ (t) )) (16)
where R θ is the Matérn correlation and θ = (φ, κ). Since the variance parameter is constant and following [START_REF] Bachoc | Parametric estimation of covariance function in Gaussian-process based Kriging models[END_REF], the optimization of the variance parameter σ 2 β can be carried out separately with the correlation parameters φ and κ. Therefore,

σ 2,(t+1) β = Tr(R -1 θ E(ββ T | y, Φ (t) )) d θ (t+1) = arg max θ ln(|R -1 θ |) -d ln(Tr(R -1 θ E(ββ T | y, Φ (t) ))). ( 17
)
The solution to the optimization problem in equation ( 17) cannot be done analytically; therefore, numerical optimization algorithms are used. This study uses the quasi-Newton method L-BFGS-B to optimize the parameters. Given the difficulties in estimating Matérn parameters [START_REF] Kaufman | The role of the range parameter for estimation and prediction in geostatistics[END_REF], we a priori fix the smoothness parameter as 3 2 , which gives the classical 3 2 -Matérn covariance function.

Conditional autoregressive model

The M-step in equation ( 10) requires the inversion of the covariance matrix, which can be challenging for large matrices. This problem is wildly discussed in Gaussian processes literature [START_REF] Ambikasaran | Fast direct methods for gaussian processes[END_REF][START_REF] Storkey | Truncated covariance matrices and toeplitz methods in gaussian processes[END_REF]. Therefore, it can be numerically advantageous to parameterize the precision matrix (inverse of the covariance matrix) instead of the covariance matrix. This is motivated by the fact that the precision matrix P θ = Σ -1 θ can be approximated by a sparse matrix [START_REF] Tajbakhsh | On the theoretical guarantees for parameter estimation of gaussian random field models: A sparse precision matrix approach[END_REF]. In fact, the off-diagonal elements of the precision matrix correspond to the conditional covariance between two variables given the remaining variables. Therefore, conditionally independent variables have zero values in the precision matrix.

Gaussian Markov random fields (GMFs) are wildly used in spatial statistics [START_REF] Cressie | Statistics for spatio-temporal data[END_REF]. GMFs models have a Markov property making them computationally and theoretically suitable [START_REF] Rue | Fast sampling of gaussian markov random fields[END_REF]. Furthermore, [START_REF] Rue | Fitting gaussian markov random fields to gaussian fields[END_REF] demonstrated that a GMF model can approximate a Gaussian field with a Matérn correlation function and other families of correlation functions. Conditional autoregressive (CAR) models are classes of GMFs with well-defined joint Gaussian distribution [START_REF] Cressie | Some diagnostics for markov random fields[END_REF]. This subsection will study cases where the coefficients β have the CAR model property. The joint distribution of a CAR is expressed as

β ∼ N (0, τ 2 (I d -αH) -1 Φ). ( 18
)
The distribution of β depends on unknown parameters α and τ 2 , and many types of CAR models depend on the choice of the matrix H and Φ. Following [START_REF] Besag | Bayesian image restoration, with two applications in spatial statistics[END_REF], in this study, we consider the Weighted CAR (WCAR) model where

Φ = diag(|N 1 | -1 , ..., |N d | -1 ) (19)
where

|N i | is the number of neighbors of location i and H = aij |Ni| d×d ; i, j = 1, ..., d
, where a ij is the (i, j) element of the adjacency matrix A = (a ij ) d×d , where a ij = a ji = 1 if and only if location i and j are neighbors and otherwise

a ij = 0. Putting P θ = τ -2 (I d -αH)Φ -1
, the second part of the M-step in the equation ( 11) becomes

θ (t+1) = arg max θ ln(|P θ |) -Tr(P θ E(ββ T | y, Φ (t) )) (20)
where θ = (τ 2 , α).

As for the Matérn covariance, the solution to the optimization problem (20) cannot be done analytically, and the numerical optimization algorithm L-BFGS-B is used. Note that the optimization of the variance parameter τ 2 can also be carried out separately with the parameter α.

Remark that this leads to a spatial extension of the fused Ridge method proposed in [START_REF] Goeman | Autocorrelated logistic ridge regression for prediction based on proteomics spectra[END_REF]. When α = 1, we obtain

1 τ 2 β T Φ -1 (I d -αH)β = 1 2τ 2 (i,j)|aij =1 (β i -β j ) 2 . ( 21
)
This shows that any spatial coefficient variations will be penalized when solving (7). In this case, replacing the L2 norm with the L1 norm leads to the fused LASSO method proposed in [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF]. However, the matrix (I p -αH) is semi-positive definite when α = 1 and thus Σ θ is degenerate. Hereafter we impose the constraints |α| < 1 to ensure that the precision matrix is positive definite. Another strategy would consist of adding a regular Ridge penalty (e.g., the discussion in van Wieringen ( 2015)).

Simulation study

In this section, a simulation study is conducted to assess the performance of the proposed method for estimating model parameters for the three cases: diagonal, Matérn, and CAR.

Setup

This study focuses on using the proposed method for spatial applications. Therefore, we consider a 15 × 15 regular spatial grid in a square domain [1, 15] 2 where each location j has a covariate x j . We generate X = (x ij ) n×d of n independent and identically distributed observations from a multivariate normal distribution with zero mean and a Matérn covariance with some arbitrary parameters (σ 2 x , φ x , κ x ) = (6, 2, 3/2). Then, the coefficients β, kept the same for all observations, are simulated using either the diagonal, Matèrn, or CAR case. Finally, for a given σ 2 , Y is simulated from the normal distribution according to equation (3).

The parameters chosen for each case are:

• Diagonal: σ 2 = 36 and σ 2 β = 7 • Matérn: σ 2 = 36, σ 2 β = 0.1 and φ = 4 • CAR: σ 2 = 36, τ 2 = 1 and α = 0.9

The parameters are chosen so that the results of the three methods are comparable. For the CAR model, we consider four neighbors to construct the adjacency matrix, and we chose α = 0.9 to sufficiently smooth the resulting coefficients.

The EM algorithm is initialized with an arbitrary set of parameters, and the E-step and M-step are repeated until no further improvement can be made to the likelihood value or to limit the computational cost until a maximum number of iterations is reached. The computation time for one iteration on an i5-7500 CPU and 16Go computer is 0.16, 3, and 1.8 seconds for diagonal, Matérn, and CAR, respectively.

Results

At first, one simulation is done for each case (diagonal, Matérn, and CAR) with n = 800. The parameters are estimated using the EM algorithm presented in the previous section. Figure 1 shows the first simulation results. Left panels correspond to the true β, and right panels correspond to the estimated β using the EM algorithm. For all the cases, the EM algorithm does well in estimating the parameters, especially the variance σ 2 .

To assess the influence of the sample size on the estimations, for each case, we perform 100 independent random simulations for each sample size varying from 50 to 850. For each simulation, the EM algorithm is used to estimate the parameters. Figure 2 shows the normalized root mean square error N RM SE β and N RM SE y for the three cases where

N RM SE β = 1 d d j (β j -βj ) 2 σβ N RM SE y = 1 n n i (y i -ŷi ) 2 σy (22)
where βj and ŷi are the estimated β j and y i and σβ and σy are the sample standard deviation of β and y, respectively. N RM SE y is calculated in a test set (which is not used in the estimation) of size n = n 2 . For the three cases, N RM SE β and N RM SE y decrease as the sample size increases.

To evaluate the parameter estimates, we compare the EM estimates with the maximum likelihood estimates of the parameters, hereafter referred to as MLE, knowing the true β. More precisely, the MLE estimates are defined as

Θ MLE = arg max Θ - 1 2 ln(|Σ θ |) + β T true Σ -1 θ β true + n ln(σ 2 ) + 1 σ 2 y -Xβ true 2 + C (23)
where β true is the true β simulated for each case with the parameters given in section 3.1. Along with the sample size, we are also interested in how the estimates behave when varying the dimension of X, d, and the variance parameter σ 2 . Note that in practice, Θ MLE cannot be found directly, given that the true β is not observed (latent variable). Therefore, we expect the EM algorithm to provide less accurate estimates than MLE. However, we expect that by varying the sample size, the dimension, and the variance σ 2 , the estimations asymptotically will be close to MLE estimates.

Figures 3, 4 and 5 show boxplots of EM (red) and MLE (blue) estimates for the diagonal, Matérn and CAR cases as a function of sample size, dimension d, and variance σ 2 . For the diagonal case, the estimate of σ 2 seems to converge to the true value of the parameter (blue line) when the sample size n increases as it does in the usual linear regression model. Note that the estimate of the spatial variance sigma 2 β does not seem to converge to the true value of the parameter as the sample size increases, but when n is large enough, EM and MLE seem to provide similar results. This is not unexpected since both methods are based on a single sample of the d-dimensional field β. As expected, the dimension d also affects the estimate of the parameter σ 2 β , which converges towards the true value as d increases; however, no significant change is observed for σ 2 when d increases. The effect of the variance σ 2 on the estimation of σ 2 β is small, and we observe that for σ 2 larger than 100, the EM and MLE tend to underestimate σ 2 β . Similar behavior can be observed for the Matérn case: the variance parameter σ 2 seems to converge towards the actual value with increasing sample size. However, there is no significant change in the other parameters (the variance σ β and the range φ). The dimension d mainly influences the parameters σ β and φ, which describe the spatial structure of the d-dimensional field β, and as d increases, the estimates converge to the actual values. As for the diagonal case, the EM algorithm underestimates the parameters σ β and φ when the variance σ 2 increases. Finally, for the CAR case, the sample size influences the parameters σ 2 and τ 2 , but only slightly the correlation parameter α, which is mainly influenced by the dimension d. The variance σ 2 has a significant influence on τ 2 , but only a small one on α. To summarize:

• The sample size n mainly influences the estimation of the variance of the residuals σ 2 The blue line corresponds to the true value of the parameter σ 2 , σ 2 β , and φ, which are equal to 36, 0.1, and 4, respectively.

σ 2 =36, σ β 2 =7 σ ^2 =36, σ ^β 2 =8
• The parameters which describe the spatial structure of β are mainly influenced by the dimension d

• As the variance σ 2 increases, EM underestimates the parameter σ 2 β of the diagonal and Matérn case, and the range parameter φ

• EM estimates are close to MLE estimates in most cases when the sample size and the dimension d are large enough and the variance σ 2 is small Another interesting aspect that needs to be studied is when the coefficients β are simulated using one covariance and estimated using another covariance model. To do that, we perform 100 independent simulations of β using the Matérn covariance function, and we estimate the parameters using the three cases: diagonal, CAR, and Matérn. Figure 6 shows the results of N RM SE β and N RM SE y of the experiment. It is clear that using the Matérn covariance for the estimation gives better results in terms of N RM SE β . Not surprisingly, the diagonal case is the worst model for estimating the coefficients. However, in terms of N RM SE y , there is a small difference between the three methods. The blue line corresponds to the true value of the parameter σ 2 , σ 2 β , and α, which are equal to 36, 1, and 0.9, respectively.

Application

The proposed method is applied to the problem of predicting the significant wave height (H s ) at a location in the Bay of Biscay using wind conditions over the North Atlantic (figure 7), where the significant wave height is the average height of the highest third of the waves, a key measure of wave height that provides information about wave energy. The data used for H s comes from the Homere hindcast database [START_REF] Boudière | A suitable metocean hindcast database for the design of marine energy converters[END_REF], and the wind data comes from Climate Forecast System Reanalysis (CFSR) [START_REF] Saha | The ncep climate forecast system reanalysis[END_REF]. The wind data are pre-processed before being used as a predictor (see [START_REF] Obakrim | Statistical modeling of the space-time relation between wind and significant wave height[END_REF] for the pre-processing procedure). We consider 23 years of H s and wind data from 1994 to 2016 with a temporal resolution of 3 hours.

The regression problem is of the form where X j (t) is the predictor at time t and location j defined as

H s (t) = d j=1 X j (t)β j + (t) t = 1, ..., n (24)
X j (t; t j , α j ) = 1 2αj +1 t-tj +αj i=t-tj -αj W 2 j (i), (25)
t j + α j + 1 ≤ t ≤ t j -α j + n
where W j is the projected wind (figure 7) defined as

W j = U j cos 1 2 (b j -θ j ) (26)
U j is the wind speed, b j is the great circle bearing, and θ j is the wind direction at location j. α j controls the length of the time window, and t j is the mean travel time of waves which are estimated using the maximum correlation between H s and the predictor ( tj , αj ) = arg max tj ,αj corr(H s , X g j (t j , α j )) .

Let X = X 1 , ..., X d be the predictor which has the size 67088 × 5651. Since the predictor has a spatial structure. It is reasonable to assume that the coefficients β also have a spatial structure so that nearby locations have close contributions to the waves at the target point. This assumption is equivalent to suppose that β ∼ N (0, Σ θ ). For the covariance Σ θ , 1: Quantitative comparison of the diagonal, Matérn, and CAR methods in the validation set using the correlation (r), root mean square error (RMSE), and bias.

we will consider the cases of Matérn and CAR. For comparison, we also consider the diagonal case even though it does not consider any structure between coefficients. The model's parameters (equation 24) are estimated using data from 1994 to 2013, and the model is evaluated in terms of correlation, RMSE, and bias, using a validation set from 2014 to 2016. Figure 8 shows the results of estimating β and the covariance parameters using the EM algorithm when the covariance structure is assumed to be diagonal, Matérn and CAR. Not surprisingly, the coefficients estimated with the diagonal covariance show no physical spatial structure. Therefore, the assumption that close locations have close coefficients cannot be taken into account using the diagonal case. This motivates using the Matérn and CAR covariances. The Matérn and CAR covariances give the smoothest coefficients with a clear spatial structure. In addition, locations close to the target point have larger coefficients. Therefore, the obtained coefficients are more physically interpretable and take into account our assumption about the covariance. Note that the CAR method is less expensive numerically than the Matérn, which involves inverting the covariance matrix at each iteration of the optimization algorithm used in the M-step.

Table 1 shows the results of the quantitative comparison between the three methods for predicting significant wave height in the validation set using correlation (r), root mean square error (RMSE), and bias. In terms of correlation and RMSE, the diagonal method is the less accurate method. Therefore, adding the spatial structure in the covariance is advantageous in predicting the significant wave height. The CAR and Matérn methods lead to close results regarding r, RMSE, and bias.

Conclusions

This study proposed an EM algorithm for estimating generalized Ridge regression with spatial covariates. We have studied three cases: the diagonal, Matérn, and the CAR case. A simulation study is carried out to evaluate the performance of the algorithms, and the EM algorithm successfully estimates the parameters in all cases. We have studied the influence of the sample size, dimension of X, and the variance σ 2 on the estimation. The sample size mainly influences the variance parameter σ 2 . The range parameter of the Matérn and correlation parameter of the CAR are mainly influenced by dimension d.

The proposed method is applied to the problem of downscaling the significant wave height in the Bay of Biscay using wind conditions over the North Atlantic. The Matérn method gives smooth coefficients with a clear spatial structure; however, the CAR method slightly outperforms the Matérn method in terms of RMSE. The Matérn covariance is clearly a better choice for spatial applications. However, estimating the parameters requires the inversion of the covariance matrix at each iteration of the optimization method in the M-step, which may be a computational bottleneck in many applications. To address this issue, instead of parameterizing the covariance matrix, one can parameterize the precision matrix directly as we did with the CAR method.

Appendices A Comparison between cross-validation and EM

As stated in section 2, the EM algorithm can be used as an alternative for cross-validation for estimating Ridge regression. In this section, we perform a simulation study to compare the two approaches and use the same simulation procedure discussed in section 3.1. Given the same covariates X (presented in section 3.1) we perform 50 independent random samples of coefficients β using the diagonal method (with parameters σ 2 = 36 and σ 2 β = 7). For each simulation, we estimate the coefficients using the EM algorithm and the cross-validation method. Figure 9 shows the box plot of N RM SE β and N RM SE y . The EM algorithm outperforms cross-validation in estimating the coefficients β and predicting y. The comparison we performed here is for the Gaussian case; therefore, it is straightforward that the EM algorithm will outperform cross-validation. To see how the two approaches behave in the non-Gaussian case, we simulate the response variable Y using the model Y = Xβ + , where ∼ U (2, 30) (A.1) Where U (2, 30) is the uniform distribution on the interval [2,30]. Figure 10 shows the estimation results using the EM algorithm and cross-validation. The EM algorithm still outperforms cross-validation in both N RM SE β and N RM SE y ; however, the difference between the two methods here is small than in the Gaussian case.

B The case where β has a non-zero mean

In this section, we consider the case where β has a non-zero mean as defined by the hierarchically model

β ∼ N (µ ξ , Σ θ ) Y | β, Θ ∼ N (Xβ, σ 2 I n ) (B.1)
where Θ = (σ 2 , µ ξ , θ).

Figure 1 :Figure 2 :Figure 3 :Figure 4 :

 1234 Figure1: Simulation results for the three cases (diagonal, CAR, and Matérn). The left panels correspond to the true β coefficients with the true parameters given in section 3.1, and the right panels correspond to the β estimated when the sample size n = 800.

Figure 5 :

 5 Figure5: Estimated parameters in the case where the covariance of β is the CAR as a function of the sample size, the dimension of X, d, and the variance σ 2 . Red boxes correspond to EM estimates and the blue ones to MLE estimates. The blue line corresponds to the true value of the parameter σ 2 , σ 2 β , and α, which are equal to 36, 1, and 0.9, respectively.

Figure 6 :

 6 Figure 6: Results of the estimations when the true beta is simulated from Matérn with the parameters σ 2 = 36, σ 2 β = 0.1 and φ = 4 and sample size n = 800. The left panel correspond to N RM SE β and the right one for N RM SE y .

Figure 7 :

 7 Figure 7: CFSR projected wind in the North Atlantic in 1994-01-01 00h:00. The black point represents the target point.

  ^2 =0.1, τ ^2 =0.00001, α ^ =0.99

Figure 8 :

 8 Figure 8: The coefficients β estimated using the EM algorithm with diagonal, Matérn, and CAR covariance.

Figure 9 :Figure 10 :

 910 Figure 9: Results of estimating Ridge regression with the EM algorithm and 10-fold cross-validation in the Gaussian case.

The complete log-likelihood is expressed as ln p(y, β; Θ) = ln p(y | β; σ 2 ) + ln p(β; θ)

Where C is a constant. In the M-step, the quantity Q(Θ|Θ (t) ) is maximized with respect to the parameters Θ.

• E-step:

The posterior distribution of the latent variable β is a normal distribution with mean µ β|y and covariance matrix Σ β|y such that

• M-step:

The maximization step computes


Said Obakrim
email: said.obakrim@univ-rennes1.fr

Valérie Monbet
email: valerie.monbet@univ-rennes1.fr

Nicolas Raillard
email: nicolas.raillard@ifremer.fr

Pierre Ailliot
email: pierre.ailliot@univ-brest.fr

LEARNING THE SPATIO-TEMPORAL RELATIONSHIP BETWEEN WIND AND SIGNIFICANT WAVE HEIGHT USING DEEP LEARNING

Keywords: Wind fields, Significant wave height, Convolutional Neural Networks, Long Short-Term Memory

Ocean wave climate has a significant impact on near-shore and off-shore human activities, and its characterisation can help in the design of ocean structures such as wave energy converters and sea dikes. Therefore, engineers need long time series of ocean wave parameters. Numerical models are a valuable source of ocean wave data; however, they are computationally expensive. Consequently, statistical and data-driven approaches have gained increasing interest in recent decades. This work investigates the spatio-temporal relationship between North Atlantic wind and significant wave height (H s ) at an off-shore location in the Bay of Biscay, using a two-stage deep learning model. The first step uses convolutional neural networks (CNNs) to extract the spatial features that contribute to H s . Then, long short-term memory (LSTM) is used to learn the long-term temporal dependencies between wind and waves.

Introduction

Characterisation of wave climate is required for many marine applications, such as the design of coastal and offshore structures and the planning of ship operations. Wind waves are generated by the surface wind, with local wind creating the wind sea and wind from distant areas creating waves that propagate and form swells ( [START_REF] Young | Wind generated ocean waves[END_REF]). Waves in the Bay of Biscay depend on both local and large-scale wind conditions in the North Atlantic [START_REF] Charles | Climate change impact on waves in the Bay of Biscay, France[END_REF]]) ; however, swells generally dominate the sea state. Swells travel large distances and take up to five days to cross the Atlantic from Cape Hatteras to the Bay of Biscay ( [START_REF] Ardhuin | Wind Waves[END_REF]]). Consequently, waves observed at a given location depend on wind conditions over the North Atlantic in a time window of several days, and it is challenging to reproduce this complex spatio-temporal relationship using machine learning. The goal of this work is to propose a deep learning approach that learns this relationship.

The advantage of deep learning methods ([Ian Goodfellow et al. (2016)]) lies in their ability to build hierarchical representations of predictors. In particular, in the case of spatial data, convolutional neural networks (CNNs) allow to learn complex spatial features from the data ([Jiuxiang Gu et al.(2018)]). Moreover, long short memory (LSTMs) ( [START_REF] Hochreiter | Long Short-Term Memory[END_REF]) have proven to be very successful in predicting time series and sequence data. In this work, we propose a non-expensive data-driven approach that learns the underlying spatio-temporal structure of the relationship between wind and waves using a two-stage model based on CNNs and LSTM. This paper is organised as follows. Section 2 presents the problem of downscaling of ocean waves and related works. Section 3 describes the data used in this work. Section 4 presents the proposed two-stage model, the architecture, and

Problem statement and related work

The problem of improving the spatial resolution of climate variables is known under the name of downscaling ( [START_REF] Maraun | [END_REF]). Downscaling approaches attempt to construct a link, either numerical or statistical, between large-scale and local-scale variables. The advantage of statistical downscaling (SD) over numerical models is primarily in terms of computational efficiency. A rigorous comparison of the two approaches can be found in ( [START_REF] Wang | Dynamical versus statistical downscaling methods for ocean wave heights[END_REF], [START_REF] Laugel | Wave climate projections along the French coastline: dynamical versus statistical downscaling methods[END_REF]).

In the case of ocean waves, wind [START_REF] Obakrim | Statistical modeling of the space-time relation between wind and significant wave height[END_REF]]) or sea level pressure (SLP) ( [START_REF] Camus | A method for finding the optimal predictor indices for local wave climate conditions[END_REF]]) are commonly used to downscale ocean wave parameters. However, in order to establish a link function between the wind (or SLP) and the local ocean wave parameters, it is necessary to consider a large spatial and temporal coverage and, consequently, a large number of potential explanatory variables that are highly correlated. Some methods determine the wave generation area for any ocean location worldwide. For example, ESTELA [START_REF] Pérez | ESTELA: a method for evaluating the source and travel time of the wave energy reaching a local area[END_REF]) is a numerical model that uses the spectral information to select the fraction of energy that travels to the target point from selected source points. The ESTELA method can be used to design statistical downscaling methods. For instance, [START_REF] Camus | A method for finding the optimal predictor indices for local wave climate conditions[END_REF]] and [START_REF] Hegermiller | A multimodal wave spectrumbased approach for statistical downscaling of local wave climate[END_REF]] used the ESTELA method to define the predictors used in their SD model. [START_REF] Obakrim | Statistical modeling of the space-time relation between wind and significant wave height[END_REF]] proposed a data-driven approach that determines the wave generation area by estimating the travel time of waves, generated in each considered sources point, that reach the target point. Then, the predictors were defined based on the wave generation area and finally a SD model based on weather types was built.

As far as we know, the existing methods for SD of ocean wave parameters define a priori the spatio-temporal structure of the predictors, and then the SD model is built using these predictors. The aim of this study is to propose a deep learning approach that automatically learns the spatio-temporal relationship between wind and waves.

Data preparation

The Climate Forecast System Reanalysis (CFSR) ( [START_REF] Saha | The NCEP climate forecast system reanalysis[END_REF]]) hourly wind data is considered in this study as a predictor. CFSR is a global reanalysis developed by the National Centers for Environmental Prediction (NCEP) that covers the period from 1979 to the present with an hourly time step and a spatial resolution of 0.5°by 0.5°. The historical H s data is extracted from the hindcast database HOMERE ( [START_REF] Boudière | Boudière, Edwige and Maisondieu, Christophe and Ardhuin, Fabrice and Accensi, Mickaël and[END_REF]]) at the target location with spatial coordinates (45.2°N, 1.6°W) located in the Bay of Biscay. The temporal resolution of both wind and H s data is up-scaled to 3-hourly data. The period from 1994 to 2016 is considered in this study, leading to a dataset with n = 67208 observations. Instead of using both zonal and meridional components as a predictor, we use the projected wind ( [START_REF] Obakrim | Statistical modeling of the space-time relation between wind and significant wave height[END_REF]]) defined, at each location j and time t, as

W j (t) = U j (t) cos 2 1 2 (b j -θ j (t)) (1)
where W j (t) is the projected wind, U j (t) the wind speed, θ j (t) the wind direction and b j is the great circle bearing from the source point j to the target point. Under the assumption that waves travel in great circle paths, grid points whose paths are blocked by land are neglected (Figure 1). Therefore, we define the global predictor at time t as

X (g) (t) = (W 2 1 (t), ..., W 2 p (t)) (2)
where p = 5651 is total number of grid points.

Following [START_REF] Obakrim | Statistical modeling of the space-time relation between wind and significant wave height[END_REF]]), in order to capture the wind sea, we also define the local predictor as

X ( ) (t) = {U (t), U 2 (t), U 3 (t), U 2 (t)F (t), U (t -1), U 2 (t -1), U 3 (t -1), U 2 (t -1)F (t -1)} (3)
where U (t) is the wind speed at the target point and F (t) is the fetch length at time t, calculated as the minimum of the distance from the target point to shore in the direction from which the wind is blowing and 500km. The fetch has an important effect on wind sea characteristics ( [START_REF] Ardhuin | Wind Waves[END_REF]]); therefore, it is commonly used to construct empirical wind wave models.

Proposed methodology

As mentioned in the last section, state of the art statistical methods for downscaling wave parameters usually use a pre-processing step to create features that take into account the wave generation area. In this study, we propose a deep learning approach that automatically extracts these features. Since waves may take several days to reach the target point, the history and current wind can be used to predict H s . An example of this type of model could have the following form

H s (t) = f (X (g) (t -t max ), ..., X (g) (t)) (4)
where, t max can be interpreted as the maximum travel time of the waves and will be referred to as such in the following. However, this approach can be computationally challenging given the dimension of the predictor (5651 in our case).

Instead, in this study we propose to use current wind conditions to estimate current and future H s .

In order to describe the complex spatio-temporal relationship between wind and H s , we propose the following two-stage model

1 st stage: [H s (t|X (g) (t)), ..., H s (t + t max |X (g) (t))] = f (X (g) (t)) + (t), f : R p → R tmax 2 nd stage: H s (t) = g(X (g) (t), f (X (g) (t -t max )), ..., f (X (g) (t))) + (t), g : R tmax * tmax+8 → R (5)
where the notation H s (t 1 |X (g) (t 2 )) represents the contribution of wind conditions at time t 2 in H s at time t 1 . and are the errors of the 1 st stage and 2 nd stage, respectively. The 1 st stage estimates the current and future H s using current wind conditions. The 2 nd stage estimates H s using the past predictions obtained from the 1 st stage. Along with the local predictor X (g) , the input for the 2 nd stage is a t max * t max matrix of the form

   Ĥs (t -t max |X (g) (t -t max )) . . . Ĥs (t|X (g) (t -t max )) . . . . . . . . . Ĥs (t|X (g) (t)) . . . Ĥs (t + t max |X (g) (t))    (6)
where Ĥs (t 1 |X (g) (t 2 )) represents the prediction, obtained from the 1 st stage, of the contribution of wind conditions at time t 2 in the H s at time t 1 . When t 1 = t 2 , this prediction represents the wind sea (first column of the matrix in equation ( 6)); for t 1 > t 2 , on the other hand, the prediction represents the H s caused by swells.

The general structure of the model is shown in Figure 2. The 1 st stage consists of a series of 3*3 convolutions followed by the ReLU activation function, 2*2 max pooling layer, Batch Normalisation, then a flatten followed by a dense layer. The 2 nd stage starts with an LSTM layer that learns the long-term dependencies of the (t -t max , ..., t) outputs of the 1 st stage. The output of the LSTM layer is then concatenated with the local predictor X l and fed into two fully

M SE(1 st stage) = 1 t max tmax i=0 1 n -t max -1 n-tmax t=1 (H s (t + i) -Ĥs (t + i|X (g) (t))) 2 M SE(2 nd stage) = 1 n n t=1 (H s (t) -Ĥs (t)) 2 (7)
Where n is the total number of observations and Ĥs is the prediction of H s . The Keras framework with Tensorflow backend ( [START_REF] Chollet | Keras[END_REF]]) is used in this work to train the model, on a Nvidia K80s GPU using the Adam optimiser ( [START_REF] Diederik | Adam: A Method for Stochastic Optimization[END_REF]]) and mini batches of 64.

Results

The period from 1994 to 2011 is used to train the two-stage model and the period from 2012 to 2014 serves as the validation period. The measures chosen in this paper to validate the analysis are the correlation coefficient (r), the root mean square error (RMSE) and the bias. Different values for the maximum travel time of waves t max are tested, and the results of k-fold cross-validation (with k = 5) are shown in Figure 3. The RMSE stabilises approximately at t max = 30 × 3h, which corresponds to about 3.3 days, and the gain is substantial compared to using t max = 5. This means that wind conditions over a time window of at least 3.3 days must be considered to characterise the wave climate at the target location. In the following, the value of t max is chosen equal to 30.

Figure 4 shows the scatter plot of observed versus predicted values of H s using the two-stage model ( 5). The RMSE in the validation period is equal to 0.21m for an H s of mean 1.9m and standard deviation 1.1m. The model performs well in predicting H s and accounts for both wind and swell. The validation measures in the calibration and validation periods are almost the same. This means that the model does not overfit the training data and generalises well the relationship between wind and waves. Furthermore, the seasonality of H s is well captured by the two-stage model, as can be seen in Figure 5.

A comparison of the two-stage model with two other statistical approaches is done in Table 1. The first approach, described in ( [START_REF] Obakrim | Statistical modeling of the space-time relation between wind and significant wave height[END_REF]]), is based on weather types ( [START_REF] Maraun | [END_REF]). As for the present work, the local and global predictors were considered. However, in order to reduce the dimension of the predictor a single predictor is extracted at each spatial location j to predict H s at time t. It is defined a priori as

X (g) j (t; t j , α j ) = 1 2αj +1 t-tj +αj i=t-tj -αj W 2 j (i), (8)
t j + α j + 1 ≤ t ≤ t j -α j + n where t j is the travel time of waves, α j controls the length of the time window, and W j is the projected wind at location j. The parameters t j and α j were estimated using the maximum correlation between h s and the global predictor. The second method ( [START_REF] Michel | Deep learning for statistical downscaling of sea states[END_REF]]) uses CNNs to predict H s using the same predictors as in [START_REF] Obakrim | Statistical modeling of the space-time relation between wind and significant wave height[END_REF]). Thus, the main difference with the approach proposed in this work is that the temporal dimension of the global predictor is reduced a priori using the preprocessing step based on maximum correlation described above. The numerical results in Table 1 indicate that the two-stage model significantly outperforms the other two methods in term of the validation measures.

Conclusion

In this study, a two-stage model based on deep learning is proposed to predict H s using wind conditions. The model is capable of learning automatically the underlying spatio-temporal structure of the relationship between wind and waves. The model does well in predicting H s and is computationally inexpensive (about 5min using a computer of 30GB RAM, 2 cores CPU, and a 16GB GPU). The proposed methodology is based on two stages which are trained separately. A natural question that arises for future work, is whether we can estimate the parameters jointly using back-propagation and eventually speed up the training process and improving the results. Future work also includes using the method to predict other sea state parameters, such as wave direction and period.

The proposed method can be used for climate and weather studies at any ocean location worldwide. For nearby locations, one can train only the 2 nd stage at each location, using the weights of one location as initialisation for the others and leaving the 1 st stage the same. The model can also learn from buoy data instead of hindcast data and eventually fill in the gaps and complete historical data.

Figure 1 :

 1 Figure 1: The projected wind defined in (2) in 1994-01-01 00h:00. The black point represents the target point

Figure 2 :

 2 Figure 2: Architecture of the two-stage model in equation (5)

Figure 3 :

 3 Figure 3: Results of cross-validation using different values of t max . The blue line represents the mean of RMSE and the red interval represents the minimum and maximum RMSE

Figure 5 :

 5 Figure 5: Time series of observed (blue line) and predicted (red line) H s in 2016

Table 1 :

 1 Comparison of the two-stage model, weather types and H-CNN methods.

		98	0.21	-0.006
	weather types	0.97	0.27	-0.03
	H-CNN	0.97	0.27	-0.04

Data Availability Statement

The hindcast data Homere is available in their website: https://marc.ifremer. fr/produits/rejeu_d_etats_de_mer_homere. The wind data is available from the CFSR website: https:// climatedataguide.ucar.edu/climate-data/climate-forecast-system-reanalysis-cfsr.

Supplementary Material For reasons of reproducibility, Python code and the processed data are available at https: //github.com/SaidObakrim/Two-stage-CNN-LSTM-.


Christian Beck
email: christian.beck@uni-tuebingen.de

Marco Grimaldo

Hender Lopez

Stefano Da Vela

Benedikt Sohmen

Fajun Zhang

Martin Oettel

Jean-Louis Barrat

Felix Roosen-Runge
email: felix.roosen-runge@mau.se

Frank Schreiber

Tilo Seydel
email: seydel@ill.eu

Immunoglobulins and Serum Albumins Consistent with a Colloid Physics Picture

Keywords: short-time protein diffusion, crowding, neutron scattering, polyclonal immunoglobulin, Bovine serum albumin, bidisperse suspensions, colloid physics, short-time diffusion, protein solution

come   L'archive ouverte pluridisciplinaire

Introduction

Understanding the diffusive transport of macromolecules in the polydisperse and crowded ensemble within the aqueous intracellular fluid of living cells is crucial to understand their function [START_REF] Maheshwari | Colloidal hydrodynamics of biological cells: A frontier spanning two fields[END_REF][START_REF] Witzel | Heterogeneities shape passive intracellular transport[END_REF] . Polydispersity [START_REF] Cho | Effect of polydispersity on diffusion in random obstacle matrices[END_REF][START_REF] Gonzalez | Impact of polydispersity and confinement on diffusion in hydrodynamically interacting colloidal suspensions[END_REF][START_REF] Ilker | Longtime diffusion and energy transfer in polydisperse mixtures of particles with different temperatures[END_REF] and crowding [START_REF] Ellis | Macromolecular Crowding: An Important but Neglected Aspect of the Intracellular[END_REF][START_REF] Höfling | Anomalous Transport in the Crowded World of Biological Cells Reports on Progress in Physics[END_REF][START_REF] Gupta | Protein Entrapment in Polymeric Mesh: Diffusion in Crowded Environment with Fast Process on Short Scales Macromolecules[END_REF][START_REF] Bucciarelli | Dramatic influence of patchy attractions on short-time protein diffusion under crowded conditions[END_REF][START_REF] Vodnala | Hard-sphere-like dynamics in highly concentrated alpha-crystallin suspensions[END_REF][START_REF] Ostrowska | Modeling crowded environment in molecular simulations Frontiers in[END_REF][START_REF] Fagerberg | Self-Diffusive Properties of the Intrinsically Disordered Protein Histatin 5 and the Impact of Crowding Thereon: A Combined Neutron Spectroscopy and Molecular Dynamics Simulation Study The[END_REF] are, thus, subject to numerous theoretical [START_REF] Makarov | Solvation and Hydration of Proteins and Nucleic Acids: A Theoretical View of Simulation and Experiment Accounts of[END_REF] , simulation [START_REF] Ostrowska | Modeling crowded environment in molecular simulations Frontiers in[END_REF][START_REF] Nawrocki | Slow-Down in Diffusion in Crowded Protein Solutions Correlates with Transient Cluster Formation[END_REF][START_REF] Feig | Crowding in Cellular Environments at an Atomistic Level from Computer Simulations[END_REF] , and experimental studies. Both living cells or even small organisms in their entire complexity [START_REF] Doster | Microscopic Diffusion and Hydrodynamic Interactions of Hemoglobin 11 in Red Blood Cells[END_REF][START_REF] Stadler | Hemoglobin Dynamics in Red Blood Cells: Correlation to Body Temperature[END_REF][START_REF] Jasnin | In Vivo Measurement of Internal and Global Macromolecular Motions in Escherichia Coli[END_REF][START_REF] Anunciado | In Vivo Protein Dynamics on the Nanometer Length Scale and Nanosecond Time Scale[END_REF][START_REF] Mamontov | Microscopic diffusion processes measured in living planarians[END_REF] as well as settings with different degrees of simplification have been explored using spectroscopic methods, such as nuclear magnetic resonance (NMR) [START_REF] Roos | NMR-detected Brownian dynamics of αB-crystallin over a wide range of concentrations[END_REF][START_REF] Roos | Coupling and decoupling of rotational and translational diffusion of proteins under crowding conditions[END_REF] , Mössbauer spectroscopy [START_REF] Nikitin | Effects of Macromolecular Crowding on Nanoparticle Diffusion: New Insights from Mossbauer Spectroscopy The Journal of[END_REF] , and fluorescence correlation spectroscopy [START_REF] Banks | Anomalous diffusion of proteins due to molecular crowding[END_REF][START_REF] Bacia | Fluorescence cross-correlation spectroscopy in living cells Nature Methods[END_REF][START_REF] Kalwarczyk | Apparent anomalous diffusion in the cytoplasm of human cells: the effect of probes' polydispersity The[END_REF] . Drastically simplified model systems with well-defined and adjustable configuration parameters can help to test models, many of which are derived from colloid physics [START_REF] Tokuyama | Dynamics of Hard-Sphere Suspensions[END_REF][START_REF] Nägele | On the Dynamics and Structure of Charge-Stabilized Suspensions[END_REF][START_REF] Banchio | Short-Time Transport Properties in Dense Suspensions: From Neutral to Charge-Stabilized Colloidal Spheres[END_REF] .

To this end, aqueous solutions consisting of a single monodisperse protein species have previously been probed, inter alia, using neutron spectroscopy and NMR [START_REF] Roos | NMR-detected Brownian dynamics of αB-crystallin over a wide range of concentrations[END_REF][START_REF] Roos | Coupling and decoupling of rotational and translational diffusion of proteins under crowding conditions[END_REF][START_REF] Le Coeur | Microscopic Protein Diffusion at High Concentration by Neutron Spin-Echo Spectroscopy[END_REF][START_REF] Roosen-Runge | Protein Self-Diffusion in Crowded Solutions[END_REF][START_REF] Rothe | Transient binding accounts for apparent violation of the generalized Stokes-Einstein relation in crowded protein solutions[END_REF] . On the diffusive short-time scale, direct interactions such as protein-protein collisions are negligible, and hydrodynamic interactions dominate which are considered crucial for biological function [START_REF] Ando | Crowding and Hydrodynamic Interactions Likely Dominate[END_REF] as well as for the quantitative understanding of long-time diffusion [START_REF] Cichocki | Long-time selfdiffusion coefficient and zero-frequency viscosity of dilute suspensions of spherical Brownian particles[END_REF][START_REF] Brady | The long-time self-diffusivity in concentrated colloidal dispersions[END_REF][START_REF] Zahn | Hydrodynamic Interactions May Enhance the Self-Diffusion of Colloidal Particles[END_REF][START_REF] Bleibel | Hydrodynamic interactions induce anomalous diffusion under partial confinement Soft Matter[END_REF] . On this short-time scale it has been found that the self-diffusion (synonymously: tracer diffusion) of well-folded proteins with a compact shape can be quantitatively understood in terms of predictions for monodisperse colloidal hard spheres [START_REF] Roosen-Runge | Protein Self-Diffusion in Crowded Solutions[END_REF] . Notably, the slowing-down of the protein self-diffusion with increasing protein volume fraction in the suspension, i.e., with increasing crowding, has been quantitatively understood for this situation. However, In actual biological systems, polydispersity is prevalent, so it is imperative to check the validity of the colloid picture in a genuinely polydisperse situation.

In this colloid picture, Stokesian dynamics simulations of binary hard spheres show that the larger component diffuses slower and the smaller component diffuses faster in a mixture than in their re-spective pure systems at the same total volume fraction [START_REF] Wang | Short-time transport properties of bidisperse suspensions and porous media: A Stokesian dynamics study[END_REF] . Realistic polydisperse simulations suggest that these rather complex systems can be mapped and understood by studying equivalent polydisperse hard-sphere models [START_REF] Ando | Crowding and Hydrodynamic Interactions Likely Dominate[END_REF] . Experimentally, the diffusion of immunoglobulin tracer proteins has recently been studied in such a naturally polydisperse and crowded setting of macromolecules [START_REF] Grimaldo | Protein short-time diffusion in a naturally crowded environment The[END_REF] and in dense phases after macroscopic phase separations [START_REF] Girelli | Molecular flexibility of antibodies preserved even in dense phase after macroscopic phase separation Molecular Pharmaceutics[END_REF] . Simulations adapted to the particular polydisperse system [START_REF] Grimaldo | Protein short-time diffusion in a naturally crowded environment The[END_REF] generalize the trend seen in binary systems [START_REF] Wang | Short-time transport properties of bidisperse suspensions and porous media: A Stokesian dynamics study[END_REF][START_REF] Grimaldo | Protein short-time diffusion in a naturally crowded environment The[END_REF] . Notably, particles smaller than the average effective radius of the macromolecular ensemble diffuse faster, and particles with a larger than average radius diffuse slower than in the monodisperse case at the same volume fraction. Computational studies showed similar effects influencing the self-diffusion close to interfaces [START_REF] Gonzalez | Impact of polydispersity and confinement on diffusion in hydrodynamically interacting colloidal suspensions[END_REF] . However, the existing explicit comparisons between experiment and colloid theory do not sufficiently address the genuine effects of polydispersity, e.g., the experimental results in Ref. 39 can still be approximately described by an effective monodisperse system due to the experimental restriction to only one type of tracer protein.

Hence, a dedicated study of a model polydisperse system with tunable tracer composition is necessary, which in its simplest case would be a bidisperse system with tunable total volume fraction and relative composition. Up to now, the separation of different diffusive contributions was possible either due to partial deuteration of the samples investigated [START_REF] Grimaldo | Protein short-time diffusion in a naturally crowded environment The[END_REF] or to the use of advanced modelling in hydrated powders [START_REF] Cisse | Dynamics of Apolipoprotein B-100 in Interaction with Detergent Probed by Incoherent Neutron Scattering The[END_REF] . Here, we test the predicted bulk behavior of an aqueous solution of two distinct tracer proteins, namely bovine serum albumin (BSA) and bovine polyclonal immunoglobulin gamma (Ig). We separate the apparent center-of-mass diffusion coefficients measured simultaneously for both proteins from the internal diffusive processes and the solvent contributions. We use both proteins in their native, protonated form and vary the protein volume fraction and ratio. This simplified setting, compared to the highly polydisperse mixtures previously studied, allows for a more quantitative comparison with results from simulations. It also circumvents the need for biological deuteration. Thus, by analyzing the data employ-ing different frameworks, we show the feasibility to investigate the short-time diffusive processes of two distinct species of label-free tracers in a solution and to separate and analyze the corresponding contributions to the scattering signal of the two proteins.

Materials and Methods

Sample preparation

Polyclonal Ig (G5009) (M Ig = 150 kDa 42 ), BSA (A3059) (M BSA = 66.4 kDa 43 ) and D 2 O were purchased from Sigma-Aldrich (now Merck KGaA) and used with no further purification. The samples were prepared by dissolving given masses of the proteins in D 2 O as described in earlier works [START_REF] Grimaldo | Diffusion and Dynamics of γ-Globulin in Crowded Aqueous Solutions[END_REF][START_REF] Beck | Nanosecond Tracer Diffusion as a Probe of the Solution Structure and Molecular Mobility of Protein Assemblies: The Case of Ovalbumin[END_REF] . The sample details can be found in Table S1.

QENS measurements

For the quasielastic neutron spectroscopy (QENS) measurements, the samples were filled in doublewalled cylindrical Aluminum sample containers with 23 mm outer diameter and a gap between the wall containing the sample fluid with a width of 0.15 mm. The containers were sealed with indium wire.

The QENS spectra were measured during the experiments 9-13-526 46 and 8-04-759 [START_REF] Grimaldo | External crowding in protein solutions[END_REF] on the neutron backscattering spectrometer IN16B [START_REF] Frick | Recent Backscattering Instrument Developments at the ILL and[END_REF] at the Institut Max von Laue -Paul Langevin (ILL) in Grenoble, France. For calibration, an empty cylinder, pure D 2 O and Vanadium were additionally measured. IN16B was used with Si(111) monochromator and analyzer crystals, corresponding to the elastic wavelength λ = 6.27 Å. The investigated q-range (0.2 Å -1 ≤ q ≤ 1.9 Å -1 ) corresponds to nanometer length scales. A phase space transformation chopper enhanced the neutron flux at the sample position at the expense of an acceptable beam divergence [START_REF] Hennig | Optimum Velocity of a Phase-Space Transformer for Cold-Neutron Backscattering[END_REF] . A standard Orange cryofurnace was employed to set the sample temperature.

Data reduction

The IN16B data were reduced and analyzed using MATLAB. The empty can contribution was sub-tracted from all samples measured. Vanadium spectra were fitted for each momentum transfer q value with a sum of two Gaussian functions to analytically account for the instrument resolution in the subsequent fits of the sample spectra. The solvent contribution was fixed based on the pure D 2 O measurements following the approach explained in Ref. 44 using the total protein volume fraction φ. Further analysis of the diffusion coefficients was performed with python3 50 using Jupyter Notebooks.

Analysis of the scattering signal with a simplified approach

For the investigated range in the energy transfer ℏω and in the momentum transfer q, the measured scattering function is dominated by the incoherent scattering signal S(q, ω) of the protein solution and can be written as the convolution of the resolution function R(q, ω) and the weighted sum of the contributions from the solvent S D 2 O (q, ω) [START_REF] Grimaldo | Diffusion and Dynamics of γ-Globulin in Crowded Aqueous Solutions[END_REF] and from the protein S Protein (q, ω) [START_REF] Grimaldo | Dynamics of proteins in solution Quart[END_REF] ,

S(q, ω) = R(q, ω)⊗[β • S Protein (q, ω) + β D 2 O • S D 2 O (q, ω)] ,
(1) where β (q) and β D 2 O (q) are scalar parameters. For globular proteins, the signal arising from the diffusive protein motions can be separated into an apparent global center-of-mass S glob (q, ω) and internal S int (q, ω) contribution:

S Protein (q, ω) = S glob (q, ω) ⊗ • • • [A 0 (q)δ (ω) + (1 -A 0 (q))S int (q, ω)] (2)
with A 0 (q) being the elastic incoherent structure factor (EISF) [START_REF] Grimaldo | High-Resolution Neutron Spectroscopy on Protein Solution Samples EPJ Web of Conferences[END_REF] . Previous studies investigated the dependence of the apparent short-time center-of-mass diffusion coefficients on the protein volume fraction for several proteins in single-component solutions [START_REF] Roosen-Runge | Protein Self-Diffusion in Crowded Solutions[END_REF][START_REF] Grimaldo | Diffusion and Dynamics of γ-Globulin in Crowded Aqueous Solutions[END_REF][START_REF] Beck | Nanosecond Tracer Diffusion as a Probe of the Solution Structure and Molecular Mobility of Protein Assemblies: The Case of Ovalbumin[END_REF] . For different globular proteins, an apparent diffusion coefficient D = D(D r , D t ), consisting of the translational D t and rotational D r diffusion, has been observed for significantly changing environments. For the time and length-scale investigated, it can be described by a Fickian diffusion process, which translates into a Lorentzian function L γ (ω) with the width γ(q) = Dq 2 . [START_REF] Girelli | Molecular flexibility of antibodies preserved even in dense phase after macroscopic phase separation Molecular Pharmaceutics[END_REF][START_REF] Beck | Nanosecond Tracer Diffusion as a Probe of the Solution Structure and Molecular Mobility of Protein Assemblies: The Case of Ovalbumin[END_REF][START_REF] Grimaldo | Dynamics of proteins in solution Quart[END_REF][START_REF] Beck | Temperature and salt controlled tuning of protein clusters[END_REF][START_REF] Sohmen | Nanosecond structural dynamics of the chaperone Hsp90[END_REF] In the energy transfer range investigated, the internal diffusive contribution of the proteins can also be described by a Lorentzian function [START_REF] Beck | Nanosecond Tracer Diffusion as a Probe of the Solution Structure and Molecular Mobility of Protein Assemblies: The Case of Ovalbumin[END_REF][START_REF] Grimaldo | Dynamics of proteins in solution Quart[END_REF] . The incoherent scattering signal can thus be approximated by

S(q, ω) = R(q, ω) ⊗ • • • β • A 0 L γ (ω) + (1 -A 0 )L γ+Γ (ω) + • • • β D 2 O • S D 2 O (q, ω)] (3)
with L γ (ω) and L Γ (ω) describing the apparent global diffusion and the internal diffusive processes, respectively. Although this approach is developed for monodisperse system, several studies have shown that it can be applied to cluster forming systems [START_REF] Beck | Nanosecond Tracer Diffusion as a Probe of the Solution Structure and Molecular Mobility of Protein Assemblies: The Case of Ovalbumin[END_REF][START_REF] Beck | Temperature and salt controlled tuning of protein clusters[END_REF] .

In this case, Equation 3is averaging on the one hand over the different global diffusive dynamics of the two proteins and on the other hand over their internal dynamics.

Separation of the BSA and Ig contributions in the scattering signal

For n different proteins in the solution, the total incoherent scattering signal from the proteins S ∑ (q, ω) can be written as a weighted sum of the different protein contributions.

S ∑ (q, ω) = n ∑ i=1 s i S i Protein (q, ω) n ∑ i=1 s i (4)
with S i Protein (q, ω) being the scattering signal of the protein i. The incoherent scattering cross-section s i of this protein is calculated as the sum s i = n p ∑ j σ j of the incoherent scattering cross-sections σ j of the atoms present in this protein multiplied by the number of this type of protein n p in solution.

To extract the individual diffusion coefficients of BSA and Ig separately from the experimental data, an advanced algorithm for fitting both the q and ℏω dependence simultaneously, have to be applied to avoid overfitting due to the spectrometer resolution, limited energy transfer, and statistical errors of the measured spectra. The fit according to Equation 3, has to be changed to

S(q, ω) = R ⊗ β (q) A 0 (q) s BSA L γ BSA (ω) + s Ig L γ Ig (ω) + (1 -A 0 (q))S int (q, ω)] + β D 2 O (q)L γ D 2 O (ω) (5)
with

γ BSA = D (BSA) exp q 2 and γ Ig = D (Ig) exp q 2 and D (Ig) exp < D (BSA)
exp . The scaling parameters s BSA and s Ig are calculated by

s Ig = σ Ig c Ig [mol/l] σ Ig c Ig [mol/l] + σ BSA c BSA [mol/l] (6)
and

s BSA = σ BSA c BSA [mol/l] σ Ig c Ig [mol/l] + σ BSA c BSA [mol/l] , ( 7
)
where σ Ig = 1011495.41 barn and σ BSA = 464377.95 barn are the incoherent scattering cross sections of Ig and BSA, calculated based on the pdb files 1IGT [START_REF] Harris | Protein Data Bank: Structure of immunoglobulin[END_REF] and 4F5S (biological assembly 1) [START_REF] Bujacz | Protein data bank: Crystal Structure of Bovine Serum Albumin[END_REF] , respectively [START_REF] Sears | Neutron Scattering Lengths and Cross Sections Neutron News[END_REF][START_REF] Jacrot | The Study of Biological Structures by Neutron Scattering From Solution Reports on Progress in Physics[END_REF] . c Ig [mol/l] and c BSA [mol/l] are the molar concentrations of Ig and BSA, respectively. The internal contribution S int (q, ω) = L (ω) described by a Lorentzian function as well as the EISF A 0 , fitted as q dependent parameter with the limitation to be monotonically decreasing, were averaged over both proteins. The very high energy resolution of IN16B is optimal to accurately determine the protein center-of-mass diffusion, which comes at the cost of a limited energy range. It has to be emphasized that this limited energy transfer available as well as the apparent background given by the solvent does not allow to separate the internal contributions of both proteins. Therefore, a simplified model is used for separating the apparent diffusion coefficients

D (BSA) exp and D (Ig)
exp of BSA and Ig, respectively, but keeping one single Lorentzian function averaging over the internal dynamics of BSA and Ig. The width of this Lorentzian function is chosen to be larger than the width of the Lorentzian function describing the apparent diffusion of BSA. Since the fit averages over the different contributions of the internal diffusion of BSA and Ig as well over the corresponding EISFs, a detailed analysis of the internal dynamics based on this data set is not possible. The goodness of fit is slightly better for the model accounting for the two distinct protein species, compared to the average population model (cf. Figure S3). This improved goodness of fit and the prior knowledge, due to the sample preparation, of the existence of the two protein species therefore justify the application of the latter model. It should be emphasized that the number of fit parameters for the bi-disperse model is increased by just one relative to the model describing only the average center-ofmass diffusion. This number of fit parameters in the fit applied to all q simultaneously is still lower than the number of fit parameters applied in the model free approach using two free Lorentzian functions for each q. By investigating the dependence of χ 2 ν as a function of the optimization parameters, the robustness of the fit can be evaluated. For the two different global fit approaches, these dependencies are shown in Figure S4.

Calculation of the theoretical diffusion coefficients

Computational studies on multi-component solutions indicate that the diffusion coefficient of a tracer changes due to the presence of a second type of tracer particle [START_REF] Wang | Short-time transport properties of bidisperse suspensions and porous media: A Stokesian dynamics study[END_REF][START_REF] Grimaldo | Protein short-time diffusion in a naturally crowded environment The[END_REF] . Therefore, to compare the experimental and simulation results for our two-component system, we use the model reported by M. Wang and J.F. Brady [START_REF] Wang | Short-time transport properties of bidisperse suspensions and porous media: A Stokesian dynamics study[END_REF] for the short-time diffusion of bidisperse suspensions of hard spheres. Specifically, we employed the reported polynomial expressions for the diffusion coefficients (for both translational and rotational). These expressions depend on the volume fraction and the composition of the mixture and are based on pairwise additive approximation (two body interactions) and semi-empirical formulas (more details can be found in the SI). For each sample condition (i.e. for a given composition and volume fraction), the translational diffusion coefficient D t and rotational diffusion coefficients D r were calculated using the volume fraction ϕ i theo = n i We use the relative deviation of the apparent center-of-mass diffusion coefficients D theo (ϕ theo ) from the monodisperse case (y = 0 and y = 1 for pure Ig and BSA, respectively) and multiplied by the corresponding experimental volume fraction dependence determined previously [START_REF] Roosen-Runge | Protein Self-Diffusion in Crowded Solutions[END_REF][START_REF] Grimaldo | Diffusion and Dynamics of γ-Globulin in Crowded Aqueous Solutions[END_REF] to obtain the theoretical apparent diffusion coefficients D theo (φ):

D (BSA) theo (φ, y) = D (BSA) theo (φ, y = 1) • D (BSA) theo (ϕ theo , y) D (BSA) theo (ϕ theo , y = 1) (8)
D (Ig) theo (φ, y) = D (Ig) theo (φ, y = 0)• D (Ig) theo (ϕ theo , y) D (Ig) theo (ϕ theo , y = 0) (9)
This calculation is necessary, since the direct conversion between the experimentally given volume fraction φ of the mixture and the effective hydrodynamic volume fraction ϕ theo determining the theoretical diffusion is only possible for the pure solutions [START_REF] Roosen-Runge | Protein Self-Diffusion in Crowded Solutions[END_REF][START_REF] Grimaldo | Diffusion and Dynamics of γ-Globulin in Crowded Aqueous Solutions[END_REF] , but not for mixtures.

By following this pathway, it is possible to express the simulation results in terms of the experimental conditions. In addition, by performing this renormalization, a possible presence of dimers of BSA [START_REF] Ameseder | Relevance of Internal Friction and Structural Constraints for the Dynamics of Denatured Bovine Serum Albumin[END_REF][START_REF] Ameseder | Homogeneous and heterogeneous dynamics in native and denatured bovine serum albumin[END_REF][START_REF] Jeffries | Preparing monodisperse macromolecular samples for successful biological small-angle X-ray and neutron-scattering experiments Nature Protocols[END_REF] and Ig 62 and slow domain motions of Ig, which are captured by the thinner Lorentzian function [START_REF] Girelli | Molecular flexibility of antibodies preserved even in dense phase after macroscopic phase separation Molecular Pharmaceutics[END_REF] , are taken into account. We note that a fraction of BSA or Ig dimers might be present in the samples, but the picture of monomers is sufficient to model the results. The effect of possible oligomers might cancel out in the comparison due to the scaling of the theory to the effective hydrodynamic size of the proteins.

Results and Discussion

Experimental average center-of-mass diffusion

To investigate the short-time self-diffusion of Ig and BSA in the same solution, we employ QENS using IN16B [START_REF] Frick | Recent Backscattering Instrument Developments at the ILL and[END_REF] , as this technique proved to be well-suited for the study of the nanosecond protein dynamics at high protein volume fractions [START_REF] Roosen-Runge | Protein Self-Diffusion in Crowded Solutions[END_REF][START_REF] Braun | Crowding-Controlled Cluster Size in Concentrated Aqueous Protein Solutions: Structure, Self-and Collective Diffusion[END_REF] (cf. Section 3.4). In a first, simpler approach, the QENS signal is analyzed for each recorded value of the momentum transfer ℏq individually. Based on Equation 3, the scattering signal of the two proteins was described by one Lorentzian function averaging over the centerof-mass diffusion and one second Lorentzian function averaging over the internal diffusive processes in the two proteins.

A representative fit result is shown in Figure 1a. The inset of Figure 1a depicts the fitted width γ versus q 2 . Although the fit averages over the two apparent diffusion coefficients associated with BSA and Ig, respectively, the relationship γ = D (av) exp q 2 (solid line) is not imposed, yet it arises naturally from the qwise fits. This indicates that on the timescale ranging from some tens of picoseconds to some nanoseconds accessible by IN16B the averaged center-of-mass diffusion of Ig and BSA undergoes a simple Fickian diffusion with an averaged apparent diffusion coefficient

D (av)
exp , comprising the translational and rotational diffusion contributions of both Ig and BSA. Due to the limited energy range |ℏω| ≤ 30 µeV of IN16B, the fitted γ(q) seemingly deviate from this relationship at the highest q (inset of Figure 1a), which can be attributed to this sampling as well as intensity limitation and is absent when spectrometers with a larger energy range are employed [START_REF] Grimaldo | Hierarchical Molecular Dynamics of Bovine Serum Albumin in Concentrated Aqueous Solution Below and Above Thermal Denaturation[END_REF][START_REF] Grimaldo | Salt-Induced Universal Slowing Down of the Short-Time Self-Diffusion of a Globular Protein in Aqueous Solution Journal of Physical[END_REF] .

The average apparent diffusion coefficients D (av) exp

obtained from samples measured at T = 295 K (green pentagons), as well as the volume fraction dependencies of the diffusion coefficients of the pure protein solutions determined earlier [BSA (red solid line) [START_REF] Roosen-Runge | Protein Self-Diffusion in Crowded Solutions[END_REF] and Ig (blue solid line) [START_REF] Grimaldo | Diffusion and Dynamics of γ-Globulin in Crowded Aqueous Solutions[END_REF] ] are shown as a function of the total volume fraction of the system φ = φ BSA + φ Ig in Figure 2b. In this case, φ Ig = c Ig • ν Ig and φ BSA = c BSA • ν BSA are the volume fractions calculated based on the partial specific volume of Ig ν Ig = 0.739 ml/g (Ref. 66) and BSA ν BSA = 0.735 ml/g (Ref. 67), respectively. Overviews of the measured samples and the corresponding parameters are given in Table S1 and Table S2. The dependence of the averaged apparent diffusion coefficient D Figure 1: Spectra with its statistical errors at q = 1 Å -1 with c BSA = 200 mg ml and c Ig = 100 mg ml at T = 295 K displayed in green. a) The total fit based on Equation 3, the averaged center-of-mass, averaged internal as well as the solvent contribution are displayed as red, blue-orange dotted, magenta and cyan lines, respectively. Inset: Width of the Lorentzian function describing the averaged center-of-mass diffusion as a function of q 2 . The solid line represents a fit of γ = D av q 2 . b) Fit result based on Equation 5 (red line), containing the contribution of the apparent center-ofmass diffusion of BSA (blue line) and of Ig (orange line), the contribution of the averaged internal diffusion (magenta line) as well as the contribution of the solvent (cyan line). The filled symbols represent the experimentally determined diffusion coefficients D of the average (green pentagons), of Ig (orange circles), and of BSA (blue squares) in the mixtures, respectively. Note that the confidence bounds on the fits depend on the mixing ratio. In samples with a very low volume fraction of one component, these result in large error bars on the symbols. ure 2b) is not monotonic, because the mixing ratio y = φ BSA /φ varies. Nevertheless, the observed average diffusion coefficients D (av) exp (φ, y) are within the limits given by D (Ig) (φ) and D (BSA) (φ) of the pure Ig and BSA solution, respectively.

Comparison with the calculated average center-of-mass diffusion

In Figure 2(a), the theoretical diffusion coefficients

D (BSA) theo and D (Ig)
theo for the individual components (BSA and Ig) in the mixture, calculated as explained in Section 3.6, are displayed for different mixing ratios y as a function of the total volume fraction φ (dash-dotted lines). The apparent diffusion coefficients for Ig in the mixture D (Ig) theo (φ) are lower than the ones for Ig in the monodisperse solution D (Ig) (φ). In contrast, for the smaller protein BSA, the apparent diffusion coefficient in the presence of Ig, D (BSA) theo (φ) is increased compared to the monodisperse case D (BSA) (φ). This observation is in agreement with several previous studies [START_REF] Wang | Short-time transport properties of bidisperse suspensions and porous media: A Stokesian dynamics study[END_REF][START_REF] Grimaldo | Protein short-time diffusion in a naturally crowded environment The[END_REF] . For the sample conditions investigated, the theoretical diffusion coefficients are displayed in addition as open symbols in Figure 2(a) (note that the mixing ratio y for each symbol is different). Based on these calculated theoretical diffusion coefficients, the averaged diffusion coefficient is calculated using the same approach as for the experimental data by approximating the sum of the two weighted Lorentzian functions with one single Lorentzian function. The calculated averaged diffusion coefficient, displayed as open green symbols in Figure 2(b), agree well with the experimentally determined ones displayed as filled green pentagons.

Individual center-of-mass diffusion of BSA and Ig

In Section 4.2, the deviation of the calculated diffusion coefficient of the protein in the mixture from the diffusion coefficient at the same volume fraction for the single component solution is described. To investigate the diffusion coefficients of BSA and Ig simultaneously, the experimental data was reanalyzed as explained in Section 3.5. Figure 1b displays an example spectrum with the fit of Equation 5extracting simultaneously the apparent global diffusion coefficients of BSA and Ig in the mixtures. The diffusion coefficients are displayed for the different samples in Figure 2(b) as filled blue squares and filled orange circles for BSA and Ig, respectively. As it can be seen, according to this analysis, the experimental diffusion coefficients in the mixtures significantly deviate from the values in the pure protein solutions at the same φ (solid lines). While the larger protein Ig is slowed down due to the presence of BSA at the same total volume fraction φ (i.e.

D (Ig) (φ) > D (Ig) exp (φ)), the smaller protein BSA is accelerated (i.e. D (BSA) (φ) < D (BSA)
exp (φ)). Even though the bi-component model seems to provide physically reasonable results, one must be aware that the data quality even under optimized experimental conditions would not allow such a separation without applying the existing knowledge on the sample composition. Overall, there is a risk of over-interpreting the data. Nevertheless, the χ 2 ν -method employed here, i.e., the maximum likelihood method for Poisson-distributed data, represents the best available statistical test to our knowledge [START_REF] Frome | Regression Analysis of Poisson-Distributed Data[END_REF] . The fact that we obtain reasonable confidence intervals stresses that we have not employed a model with redundant fit parameters.

To investigate the dependence on the mixing ratio y, the relative change of the diffusion coefficients

D = D i exp/theo -D i D i ( 10
)
for i representing BSA and Ig, respectively, is shown for both proteins in Figure 3 as a function of φ for both proteins for different y.

A systematic trend is seen for both proteins within the investigated range of y. Here, the bigger protein Ig is slowed down while the smaller protein present in solution is accelerated. It should be mentioned, that approaching the limits y = 0 and y = 1, the diffusion coefficients of the proteins in the mixture have to approach the values for the respective pure solutions:

D (BSA) exp (φ, y → 1) ! → D (BSA) (φ) (11)
D (Ig) exp (φ, y → 0) ! → D (Ig) (φ). (12)
The good agreement between the experimental D i exp and the theoretically calculated D i theo (Section 3.6) is confirmed in Figure 4, where the individual experimental diffusion coefficients for each component D i exp are plotted as a function of the calculated diffusion coefficient D i theo for each sample (symbols). We note that Ig is not spherical which might cause a small systematic error in the calculation of the apparent diffusion coefficient D i theo based on the translational and rotational diffusion coefficient. The bisector marked by the dotted line would correspond to a perfect agreement of the result from the fit of the measured spectra in terms of the bidisperse model (Equation 5) with the result from the calculation according to Wang and Brady [START_REF] Wang | Short-time transport properties of bidisperse suspensions and porous media: A Stokesian dynamics study[END_REF] . This figure also illustrates the distinct diffusion coefficients of the two components, namely Ig and BSA, in the mixture (orange circles and blue squares, respectively). We Filled blue squares and filled orange points represent the experimental diffusion coefficients of BSA and Ig, respectively. Note that the confidence bounds on the fits depend on the mixing ratio. In samples with a very low volume fraction of one component, these result in large error bars on the symbols.

have not considered the non-isotropic shape of the Ig proteins and their possible patchy interactions, although both anisotropy and charge-mediated interactions have been shown to influence protein-protein interactions in the short-time limit [START_REF] Bucciarelli | Dramatic influence of patchy attractions on short-time protein diffusion under crowded conditions[END_REF][START_REF] Skar-Gislinge | A Colloid Approach to Self-Assembling Antibodies Molecular Pharmaceutics[END_REF] . We attribute the high consistency of our results with a picture of purely hydrodynamic interactions to the fact that our Ig is polyclonal, thus featuring many different charge patterns and resulting in no net large effect of the overall patchy interactions, and to self-buffering of the protein solutions at high concentrations [START_REF] Da | Effective interactions and colloidal stability of bovine γ-globulin in solution[END_REF] .

Due to the restricted energy range of the spectrometer, a further separation of the internal dynamics contributions, which are significantly faster than the center-of-mass diffusion, is not possible. In Figure 5 the averaged EISF is shown for the different samples investigated. We used previously established models to describe the q dependence 44 :

EISF(q) = p + (1 -p) • (ϕ • A 3 j (q) + (1 -ϕ ) • A s (q)) (13) A 3 j (q) = 1 3 (1 + 2 j 0 (qa)) (14)
A s (q) = 3 j 1 (qR) qR 2

with a = 1.715 Å being the averaged distance between the hydrogens in the methyl group to fit the EISF. The dependence of p, ϕ and R on the mixing ratio y are shown in Figure 5. The width of the Lorentzian describing the averaged internal diffusion is approximated by a jump diffusion model to describe the q dependencies. The corresponding fits and their results are shown in Figure S5 in the Supporting Information. While the parameters for the EISF p and ϕ show trends as a function of y, the jump diffusion parameters, namely the residence time τ as well as the diffusion coefficient of D int , stay nearly constant at τ ≈ 0.1 ns and D int ≈ 100 Å 2 /ns, respectively, as shown in Figure S5. Also, R remains nearly constant close to the value of 10 Å. Therefore, average protein dynamics between BSA and Ig as seen by QENS appear to be rather similar. Otherwise, one would see random jumps in the internal dynamics when the mixing ratio is varied. From a biological point of view, the observation presented in this work is particularly interesting for cases in which reactions are limited or enhanced by diffusion and crowded settings [START_REF] Dorsaz | Diffusion-limited reactions in crowded environments[END_REF] . Examples are red blood cells, where the slow long-time diffusion of concentrated hemoglobin was recently suggested to be crucial in maximizing the oxygen capture at the cell level during the time spent near the alveolar sac [START_REF] Longeville | Hemoglobin diffusion and the dynamics of oxygen capture by red blood cells[END_REF]

Conclusions

Using quasi-elastic neutron scattering, we have probed the self-diffusion of proteins in crowded bidisperse supsensions at mixing compositions which were precisely set by the design of the experiment. This design constitutes a minimal model system to study effects of polydispersity in a controllable manner. We have provided a benchmark for the analysis of high resolution quasi-elastic neutron scattering spectra by accounting for two distinct proteins in aqueous solution in our model scattering function. We successfully obtained the apparent global center-of-mass diffusion coefficients for both proteins. In addition, the averaged internal diffusion contributions were separated from the signal. Within the experimental accuracy, we find quantitative agreement between the corresponding colloid model system of bidisperse hard spheres and the experimental results for the protein center-of-mass short-time self-diffusion. In particular, there are significant deviations of the diffusion for each component in the mixture when compared to the monodisperse case both in experiment and colloid theory. This short-time diffusion on the nanosecond time scale, where hydrodynamic interactions dominate, constitutes an important quantity for both the cellular function and for the calibration of long-time diffusion. Our results illustrate the predictive power of colloid hard sphere models for protein diffusion on the observation scale of short-time diffusion for the situation of two distinct protein sizes. These results contribute to a better understanding of the role of macromolecular polydispersity in living systems. This polydispersity and resulting dispersion of diffusion rates influences the diffusive transport in living cells and, thus, their function.

4 3 π R i H 3

 43 obtained employing the number density n i = c p [mg/ml]/M i W • N A , the molecular weight M i w , the hydrodynamic radius R i H for protein i and the Avogadro constant N A . The apparent center-of-mass diffusion coefficient D is calculated subsequently using the implicit relation for D = D(D r , D t )[START_REF] Roosen-Runge | Protein Self-Diffusion in Crowded Solutions[END_REF] (see SI for details).

  (av) exp on φ (green pentagon symbols in Fig-

Figure 2 :

 2 Figure 2: Diffusion coefficients as a function of the total volume fraction φ. The red and blue solid lines indicate the parametrization of the experimentally determined φ-dependence of the short time self diffusion coefficient of pure BSA and Ig, respectively. a) The dotted and dashed-dotted lines represent the short time self-diffusion of BSA and Ig, respectively, calculated as explained in the main text based on Wang and Brady 38 in mixtures containing both proteins with the mixing ratio y color-coded in the colorbar on the right. The calculated values for the samples investigated are additionally represented by open orange circles and open blue squares for pure Ig and pure BSA, respectively. b) The averaged calculated diffusion coefficient in the mixture is displayed as open green pentagrams for the sample conditions investigated. The filled symbols represent the experimentally determined diffusion coefficients D of the average (green pentagons), of Ig (orange circles), and of BSA (blue squares) in the mixtures, respectively. Note that the confidence bounds on the fits depend on the mixing ratio. In samples with a very low volume fraction of one component, these result in large error bars on the symbols.

Figure 3 :

 3 Figure 3: Relative change of the diffusion coefficients D as a function of the volume fraction φ. The corresponding mixing ratio y is color-coded. Dashed dotted lines represent the deviations predicted by Brady and Wang. The corresponding theoretical values for the specific sample conditions are additionally displayed as open blue squares and open orange circles for BSA and Ig, respectively. Blue and orange filled symbols represent the experimentally determined percentage deviations of BSA and Ig, respectively, in the mixture compared to the monodisperse solution. Note that the confidence bounds on the fits depend on the mixing ratio. In samples with a very low concentration of one component, these result in large error bars on the symbols.

Figure 4 :

 4 Figure 4: Experimentally determined diffusion coefficients D i exp as a function of the corresponding calculated diffusion coefficients D itheo according to the model by Wang and Brady[START_REF] Wang | Short-time transport properties of bidisperse suspensions and porous media: A Stokesian dynamics study[END_REF] as explained in the text. Filled blue squares and filled orange points represent the experimental diffusion coefficients of BSA and Ig, respectively. Note that the confidence bounds on the fits depend on the mixing ratio. In samples with a very low volume fraction of one component, these result in large error bars on the symbols.

  or the DNA replication, where nucleotides are diffusing

Figure 5 :

 5 Figure 5: Averaged EISF as a function of q of the different samples investigated. Fits of Equation 13 are shown as solid lines. Figures b-d represent the fit parameters p, ϕ and R as a function of the mixing ratio y. Red dashed lines are guides to the eye.

Acknowledgment

We acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG) and Agence Nationale de la Recherche (ANR-16-CE92-0009-01, Im-munoglobulinCrowding), the German Ministry of Research BMBF (05K19VTB), as well as support by the Partnership for Soft Condensed Matter (PSCM) in Grenoble. We thank M.K. Feustel née M.K. Braun for help during the experiments and the machine shop of the IAP, University of Tübingen, for producing the sample containers. We are grateful to the allocation of neutron beam time by the ILL.

Associated Content

Supporting Information Available: The Supporting Informations contain details on the calculation of the theoretical translational and rotational diffusion coefficient, as well as on the calculation of the apparent diffusion coefficients and an overview of the measured samples.

TOC Graphic