domain: sdv math stat

text:
Marco Guazzini

Alexander G Reisach

Sebastian Weichwald

Christof Seiler

spillR: Spillover Compensation in Mass Cytometry Data

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

Introduction

Mass cytometry makes it possible to count a large number of proteins simultaneously on individual cells [START_REF] Bandura | Mass cytometry: Technique for real time single cell multitarget immunoassay based on inductively coupled plasma time-of-flight mass spectrometry[END_REF][START_REF] Bendall | Single-cell mass cytometry of differential immune and drug responses across a human hematopoietic continuum[END_REF]. Although mass cytometry has less spillover-measurements from one channel overlap less with those of another-than flow cytometry [START_REF] Bagwell | Fluorescence spectral overlap compensation for any number of flow cytometry parameters[END_REF][START_REF] Novo | Generalized unmixing model for multispectral flow cytometry utilizing nonsquare compensation matrices[END_REF], spillover is still a problem and affects downstream analyses such as differential testing [START_REF] Weber | diffcyt: Differential discovery in high-dimensional cytometry via high-resolution clustering[END_REF][START_REF] Seiler | CytoGLMM: Conditional differential analysis for flow and mass cytometry experiments[END_REF] or dimensionality reduction [START_REF] Mccarthy | Scater: Pre-processing, quality control, normalization and visualization of single-cell RNA-seq data in R[END_REF]. Reducing spillover by careful design of experiment is possible [START_REF] Takahashi | Mass cytometry panel optimization through the designed distribution of signal interference[END_REF], but a purely experimental approach may be neither sufficient nor efficient [START_REF] Lun | Influence of node abundance on signaling network state and dynamics analyzed by mass cytometry[END_REF]. [START_REF] Chevrier | Compensation of signal spillover in suspension and imaging mass cytometry[END_REF] propose a method for addressing spillover by conducting an experiment on beads. This experiment measures spillover by staining each bead with a single type of antibody. The slope of the regression line between target antibody and non-target antibodies represents the spillover proportion between channels. [START_REF] Miao | Ab initio spillover compensation in mass cytometry data[END_REF] attempt to solve spillover by fitting a mixture model. Our contribution combines the solutions of [START_REF] Chevrier | Compensation of signal spillover in suspension and imaging mass cytometry[END_REF] and [START_REF] Miao | Ab initio spillover compensation in mass cytometry data[END_REF]. We still require a bead experiment, as in [START_REF] Chevrier | Compensation of signal spillover in suspension and imaging mass cytometry[END_REF], but estimate spillover leveraging a statistical model, as in [START_REF] Miao | Ab initio spillover compensation in mass cytometry data[END_REF]. Both previous versions rely on an estimate for the spillover matrix. The spillover matrix encodes the pairwise spillover proportion between channels. We avoid estimating a spillover matrix and instead model spillover by fitting a mixture model to the observed counts. Our main new assumption is that the spillover distribution-not just the spillover proportion-from the bead experiment carries over to the biological experiment. In other words, we transfer the spillover distribution to the real experiment instead of just the spillover proportion encoded in the spillover matrix. The main difference between our method and the one by [START_REF] Chevrier | Compensation of signal spillover in suspension and imaging mass cytometry[END_REF] is the average count after correcting for spillover. Our method will increases the average count after correction. By contrast, the method by [START_REF] Chevrier | Compensation of signal spillover in suspension and imaging mass cytometry[END_REF] will shrink all the counts towards zero. As a consequence, the average count will be lower after correction.

In Section 2, we present our mixture model and link it to calculating spillover probabilities for specific count values. Our estimation procedure is based on an EM algorithm and logistic regression, and implemented in our new R package spillR1 . In Section 3, we conduct experiments on simulated and real data obtained from the CATALYST R package [START_REF] Chevrier | Compensation of signal spillover in suspension and imaging mass cytometry[END_REF]. Section 4 discusses our synthetic experiments and relates our findings to CATALYST.

Methods

Example

Figure 1 illustrates our procedure using a dataset from the CATALYST package as an example. There are four markers, HLA-DR (Yb171Di), HLA-ABC (Yb172Di), CD8 (Yb174Di), and CD45 (Yb176Di), that spill over into the target marker, CD3 (Yb173Di). The markers have two names: the first name is the protein name and the second name in brackets is the conjugated metal. There are bead experiments for each of the spillover markers.

Panel A depicts the marker distributions from the beads experiment. We see that for this marker the bead experiments are high-quality as the target marker Yb173Di is concentrated around six, similarly to the experiment with real cells. This suggests that the spillover marker values can be transferred to the real experiments. Marker Yb172Di shows large spillover into Yb173Di, and suggests that the left tail of the first mode of the distribution may be attributed to that marker. The other spillover markers have low counts, making it justifiable to set some or all the low counts to zero.

Panel B has a solid black and gray curves representing our spillover probability estimates. With the smoothing parameter k = 11, we can see that the probability of spillover goes up to around 0.5. In that case, our correction step assigns around 50% of cells to spillover, and keeps the other 50% at the current value. We add a black dashed line to all plots at the position 2.725 to illustration this point. Low counts have spillover probability of one, which means that our procedure assigns them to spillover and masks them from the sample by setting them to NA values. Counts above four stem from spillover with probability zero (and from the actual target with probability one), which means that our procedure keeps them at their raw uncorrected value. The estimated curve with smoothing parameter k = 3 is more irregular, suggesting that these fluctuation may be driven by noise. Users can control the smoothing parameter k to choose the desired bias-variance tradeoff.

Panel C displays the distribution of our target marker, CD3 (Yb173Di), before and after spillover correction.

Our compensation method, spillR, masks most markers in the low counts up to two, about 25% to 50% in the medium counts range between two and four, and keeps counts above four. Larger counts are not affected by the correction, as can be seen by the overlap of both curves. In contrast, the other method, CATALYST, shifts large counts to the left, and shifts medium counts to low counts or zero counts. Zero counts (or equivalently NA counts for spillR) and mean counts are shown in Table 1.

C

Figure 1: A: Gaussian kernel density plot of target and spillover markers. B: The solid black curve represents the spillover probability estimated with our method with smoothing parameter k = 11. The grey curve is the estimate with k = 3. Larger values yield smoother spillover probability curves. C: Frequency polygons of marker CD3 (Yb173Di) comparing no correction, our method, and another method. Zero counts are not plotted. See Table 1 for the zero counts and means. Counts are arcsinh transformed with cofactor of five [START_REF] Bendall | Single-cell mass cytometry of differential immune and drug responses across a human hematopoietic continuum[END_REF]. The vertical dashed line helps to interpret the spillover correction when CD3 (Yb173Di) is 2.725. This is an interesting point as it balances between spillover counts from Yb172Di measured in beads (panel A) and the counts measured on real cells (panel C). The spillover probability is around 0.5. This means we will set about half of the counts to NA's at this point. After correction the read curve (panel C) will be adjusted downwards to the green curve (panel C) with our method spillR.

Definition of Spillover Probability and Assumptions

We observe a count Y i of a target marker in cell i. We model the observed Y i as a finite mixture [START_REF] Mclachlan | Finite mixture models[END_REF] of unobserved true marker counts Y i | Z i = 1 and spillover marker counts

Y i | Z i = 2, . . . , Y i | Z i = K with mixing probabilities π k = P (Z i = k) for k = 1, . . . , K, P (Y i = y) = K k=1 π k P (Y i = y | Z i = k).
The first mixing probability is the proportion of true signal in the observed counts. The other K -1 mixing probabilities are the proportions of spillover. The total sum of mixing probabilities equals one, k π k = 1. The total number of markers in mass cytometry panels is between 30 and 40 [START_REF] Bendall | Single-cell mass cytometry of differential immune and drug responses across a human hematopoietic continuum[END_REF], but only a small subset of three to four markers spill over into the target marker [START_REF] Chevrier | Compensation of signal spillover in suspension and imaging mass cytometry[END_REF]. So, typically

K = 1 + 3 or K = 1 + 4.
Experimentally, we only measure a sample from the distribution of Y i . The probabilities π k and true distributions P (Y i = y | Z i = k) are unobserved, and we need to estimate them from data. In many applications, the mixture components are in a parametric family, for example, the negative binomial distribution. As spillover correction is a pre-processing step followed by downstream analyses, choosing the wrong model can introduce biases in the next analysis step. To mitigate such biases, we propose to fit nonparametric mixture components. We make two assumptions that render the components and mixture probabilities identifiable:

• (A1) Spillover distributions are the same in bead and real experiments.

The distribution of Y i | Z i = k for all k > 1 is the same in beads and real cells. This assumption allows us to learn the spillover distributions of Y i | Z i = k for all k > 1 from experiments with beads, and transfer them to the experiment with real cells. This assumption relies on high-quality single stained bead experiments that measure spillover in the same range as the target biological experiment. In other words, a high-quality bead experiment for our method works best if the distribution of bead cells is similar to the distribution of real cells.

• (A2) For each cell i, the observed count Y i can only be due to one distribution.

This assumption is already implied by the statement of the mixture model. It allows us to calculate the spillover probability for a given count Y i = y from the posterior probability that it arises through spillover from markers k > 1,

P (spillover | Y i = y) = P (Z i > 1 | Y i = y) = 1 -P (Z i = 1 | Y i = y) = 1 - π 1 P (Y i = y | Z i = 1) P (Y i = y) .
To parse this calculation, recall that in mixture models the π 1 is the prior probability,

P (Y i = y | Z i = 1)
is the conditional probability given the mixture component, and the denominator P (Y i = y) is the marginal distribution. Applying Bayes rule leads to the posterior probability.

Estimation of Spillover Probability

We propose a two step procedure for estimating the spillover probability. In step 1, we estimate mixture components and mixture probabilities. We refine these estimates using the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. In step 2, we use these probability estimates to assign counts to spillover or signal.

We denote the n × K count matrix as Y = (y ik ) with real cells in the first column and beads in columns two and higher. To simplify mathematical notation but without loss of generality, we assume that the number of cells from real and bead experiments have the same n. In practice, the number of cells from bead experiments is much smaller than from real experiments. The kth column of Y contains marker counts for the kth spillover marker, which represents the empirical spillover distribution of marker k into the target marker, that is, the marker in the first column of Y.

We use the empirical weighted cumulative distribution function (CDF) with non-negative weights w i at each data point to model the marker distributions F for a marker k,

F k (y) = n i=1 w i 1 (y ik ≤ y) n i=1 w i ,
where 1(•) is the indicator function.

EM Algorithm

• Initialization: For the mixture probability vector, we assign probability 0.9 to the the target marker and divide the probability 0.1 among the spillover markers, π1 = 0.9 and πi = 0.1/(K -1) for all i > 1, and evaluate the kth mixture component using its empirical CDF with equal weights for all data points,

P (Y i = y | Z i = k) = F k (y) -F k (y -1) with w i = 1/n.
We smooth the probability mass function (PMF), P (Y i = y | Z i = k), using running medians with a fixed window size as implemented in R function runmed. The procedure is not sensitive to the choice of the initial mixture probability vector. Other settings are possible, e.g., setting all probabilities to the same value.

• E-step: We evaluate the posterior probability of a count y belonging to component k (that is, originating from marker k),

P (Z i = k | Y i = y) = πk P (Y i = y | Z i = k) K k =1 πk P (Y i = y | Z i = k )
.

• M-step: We estimate the new mixture probability vector from posterior probabilities, πk = 1 n

n i=1 P (Z i = k | Y i = y) ,
and estimate the new target marker distribution using its empirical CDF with weights set to posterior probabilities,

P (Y i = y | Z i = 1) = F 1 (y) -F 1 (y -1) with w i = P (Z i = 1 | Y i = y) .
As before, we smooth the PMF, P (Y i = y | Z i = 1), using running medians with the same window size. We keep the bead distributions, P (Y i = y | Z i = k) for all k > 1, fixed at their initial value.

To refine our estimates, we iterate over the E and M-steps until estimates stabilize. We stop iterating when π1 changes less than 10 -5 from the previous iteration. The final output is the spillover probability curve with estimates at discrete points in the support of Y i ,

P (spillover | Y i = y) = 1 -P (Z i = 1 | Y i = y).
We rely on assumption (A1) to justify updating only the distribution of the target marker. We rely on assumption (A2) to justify calculating the spillover probability from the mixture model. We refer to Appendix A for a step-by-step example of our EM algorithm.

Spillover Decision

To perform the spillover compensation, we draw from a Bernoulli distribution with the spillover probability as parameter to decide whether or not to assign a given count to spillover. We consider spillover counts as having no clear biological interpretation and mask them from our dataset while keeping all other counts. In our implementation, we choose to set spillover counts to NA instead of zero to avoid zero-inflated distributions.

Results

We first evaluate our new method spillR on simulated datasets. We probe our method to experimentally find its shortcomings. Then, we compare spillR to the non-negative least squares method implemented in the R package CATALYST on real data from the same package. All experiments and plots can be reproduced by compiling the R markdown file spillR_paper.Rmd2 .

Simulated Data

We choose three different experiments to test spillR against different bead and real cell distributions. We explore a wide range of possible parameter settings. Figure 2 has three panels, each representing one experimental setup. The first two panels test our assumptions (A1) and (A2). The third panel tests sensitivity of spillR to bimodal bead distributions. For all three experiments, we model counts using a Poisson distribution with parameter λ. We simulate 10,000 real cells with λ = 100, and 1,000 beads with λ = 70, and spillover probability of 0.1. Beads are an independent copy of the true spillover. The other parameters and statistical dependencies are specific to each experiment. The details of the generative models are given in Appendix B. We repeat each simulation 20 times and report averages over the 20 replications. Each panel of Figure 2 has two rows of plots. The plot on the first row represents the summary of the means for each experimental setup as a function of their respective parameter τ . This parameter has a different meaning in each setup. To visualize the different experiments, we summarize the full distributions with the true simulated signal mean (black), the uncorrected mean (orange), and the spillR corrected mean (green). Plots on the second row illustrate the simulated data distributions for three selected parameters τ picked from the experimental setup. The yellow density curve is the observed count Y . The black density curve is the target cell count. The blue density curve is the spillover distribution. The goal of the experiment is to estimate the mean of the black density as accurately as possible from the yellow density curve, which represent the data Y that we would observe in practice. We simulate this data ourselves with the models in Appendix B.

In the first experiment (panel A), we shift the measured beads spillover away from the true bead spillover to probe (A1). We test a wide range of bead shifts from no shift at τ = 0 to τ = -30. At τ = -30, the measured spillover (the first mode of the yellow density) is shifted away from the actual spillover (the blue density). Such low-quality beads cause both the observed and compensated mean to be below the true mean. As the beads quality improves, the compensated signal move closer to the true mean. As we increase τ the first mode of the yellow density moves towards the blue density. This represents high-quality bead experiments. In all cases, even for bad beads experiments, our compensation improves the means. Our compensation also increases means, as it should, in contrast to e.g. [START_REF] Chevrier | Compensation of signal spillover in suspension and imaging mass cytometry[END_REF].

In the second experiment (panel B), we mixed target and spillover to explore the robustness of our method with respect to our second assumption (A2). One way to think about this is that the mixture is a form of model misspecification. Our mixture model is undercomplete, which means that there are more true mixture components than we observe in the beads experiment. If τ = 0, then assumption (A2) is correct, but for τ = 0.5 the assumption (A2) is maximally violated. The true mean decrease with increasing τ . Our compensation is closer to the true mean across the tested range. At τ = 0.5 all three distributions and their means are the same.

In the third experiment (panel C), we model spillover with a bimodal distribution. Here τ is the mixing probability of the two modes. The locations of the two spillover modes are fixed. If τ = 0 or τ = 1, then spillover is unimodal. If τ = 0.5, the first mode of the bimodal beads distribution is left to the signal mode, and the second mode is to the right. The corrected mean is closer to the true mean than the uncorrected mean across the test range.

Real Data

We compare our method to CATALYST on one of the example datasets in the CATALYST package. The dataset has an experiment with real cells and a corresponding bead experiment. The experiment on real cells has [START_REF] Bendall | Single-cell mass cytometry of differential immune and drug responses across a human hematopoietic continuum[END_REF]. 5,000 peripheral blood mononuclear cells from healthy donors measured on 39 channels. The experiment on beads has 10,000 cells measured on 36 channels. They have single stained bead experiments. The number of beads per metal label range from 112 to 241.

We compare the two methods on the same markers as in the original CATALYST paper [START_REF] Chevrier | Compensation of signal spillover in suspension and imaging mass cytometry[END_REF] in their Figure 3B. In the original experiment, they conjugated three proteins-CD3, CD8, and HLA-DR-with two different metal labels. They conjugated CD8 (first row in Figure 3) with Yb174Di (Yb is the metal and 174 is the number of neutrons of the isotope) and La139Di, and similarly for the other rows. On the horizontal axis, we plot the same markers as in the original paper, CD3 and HLA-ABC. We visualize the joint distributions using two-dimensional histograms.

In all six panels (A-F), we observe that spillR compensates most strongly in the low counts. In panel C, CD3 (Yb173Di) against HLA-ABC (Yb172Di), CATALYST can be seen to compensate strongly in the middle range. It removes the spherical pattern that shows correlation between the two markers. spillR preserves this correlation structure and only masks out the lower counts of CD3 (Yb173Di). This highlights a key difference between spillR and CATALYST: spillR does not correct all counts by shrinking them, but rather removes some counts, following the idea that the distribution of the remaining counts is close to the true distribution. CATALYST follows another strategy by shrinking counts across the entire range.

The color code of the two-dimensional histograms indicates the absolute number of cells that fall into one hexagon bin. The uncorrected and spillR corrected histograms can contain different absolute number of cells because of how spillR rounds counts to integers. The raw mass cytometry data is often not count data due to proprietary post-processing of the manufacturer of the mass cytometer. That is why we convert mass cytometry data to count data, we convert the raw values to the next lower integer. The uncorrected counts do not undergo this pre-processing step. CATALYST does not perform this pre-processing step. This also explains the different patterns in panel B. spillR has horizontal strips that correspond to non-integer values not in the support of the distribution for spillR.

Discussion

The experiment for (A1) shows that the mean count after spillR correction is closer to the true mean over a wide range of bead shifts. This indicates that our method can perform well even if the bead experiments are imperfect. If the difference between distributions of beads and real cells is large, then one option is to rerun the bead experiments to reduce this gap. The experiment for (A2) shows that our method is also robust to model misspecification. Additionally, misspecification can be addressed by adding all channels if necessary. The increase in computational cost when adding channels is relatively minor as our method scales linearly in the number of spillover markers. The experiment on bimodal bead distributions shows that the mean count after correction is still closer to the true mean even with bimodal bead distributions, and also if the spillover is actually larger than the true signal.

In our comparison to CATALYST on real data, we observe the effect of the two different correction strategies. CATALYST essentially shrinks counts towards zero by minimizing a non-negative least squares objective. It assumes that spillover is linear up to counts of 5,000. The applied shrinkage is the same for low counts (e.g., below 10) and high counts (e.g., more than 100). By contrast, spillR does not require linearity of the spillover, but assumes that the distribution on the beads experiment carries over to the real cells experiment. In other words, the optimal bead experiments has the same peaks than in the real experiments.

If counts are in the spillover range (which mostly applies to low counts), they are corrected strongly and set to NA values. If counts are not in the spillover range, then they are left unchanged. Despite setting values to NA, correlations between markers are preserved. The marker correlation between HLA-ABC (Yb172Di) and CD3 (Yb173Di) illustrates this point. CATALYST removes the positive correlation, whereas spillR keeps the correlation for the higher counts. Compensation methods should try to remove spillover while keeping biological meaningful signal for unbiased downstream analyses. Further experiments on the correlation structure between these markers are necessary to resolve the discrepancy between the two methods. This is an important point as discovering correlations between markers can lead to the discovery of new clusters or signaling networks.

Another advantage of our method is the diagnostic plot of the spillover probability curve. We can judge if the curve makes sense by comparing it to the observed count and bead distributions. Methods based on non-negative least squares are harder to diagnose as they minimize a cost function with no clear biological interpretation.

Currently, we do not take advantage the target bead distribution in our estimation procedure. We only use spillover bead distributions. In future work, we aim to investigate ways to incorporate the target distribution into our estimation in a nonparametric manner. In our view, one of the biggest strengths of our current method is that it does not assume a specific parametric model for count data. We believe that this is crucial because spillover is just one step that precedes many downstream analysis steps, and avoiding the introduction of bias is thus our top priority.

A EM Algorithm Example

Here we illustrate the procedure using a numerical example that includes one target and one spillover marker. We have one data matrix Y that contains real cell counts recorded for marker 1 (column 1) and the bead counts for marker 1 when the true marker was marker 2 (column 2). In practice, Y is usually a matrix with more than two columns representing multiple spillover markers. The index i is a specific cell in beads and real cells experiment, respectively. Let's assume the following counts,

Y = (y ij ) =         3 2 5 3 17 2 3 17 2         . target <-c(3, 5, 17, 3, 17, 2) spillover <-c(2, 3, 2, NA, NA, NA) Y = dplyr::bind_cols(target = target, spillover = spillover) Y ## # A tibble: 6 x 2 ## target spillover ## <dbl> <dbl> ## 1 3 2 ## 2 5 3 ## 3 17 2 ## 4 3 NA ## 5 17 NA ## 6 2 NA
• Initialization: We initialize our EM algorithm by estimating the conditional probability of observing y given that it belongs to the target marker, and another conditional probability given that it belongs to the spillover marker,

P (Y i = 2 | Z i = 1) = 1/6 P (Y i = 2 | Z i = 2) = 2/3 P (Y i = 3 | Z i = 1) = 2/6 P (Y i = 3 | Z i = 2) = 1/3 P (Y i = 5 | Z i = 1) = 1/6 P (Y i = 17 | Z i = 1) = 2/6. n1 <-sum(!is.na(Y$target)) n2 <-sum(!is.na(Y$spillover)) F1 <-spatstat.geom::ewcdf(Y$target, weights = rep(1/n1, n1
)) F2 <-spatstat.geom::ewcdf(Y$spillover, weights = rep(1/n2, n2)) par(mfrow = c(1, 2)) plot(F1, main = "CDF of Target", xlab = "y", ylab = "F1(y)") plot(F2, main = "CDF of Spillover", xlab = "y", ylab = "F2(y)") ## target P_Y1 P_Y2 ## 1 2 0.6923077 0.3076923 ## 2 3 0.9000000 0.1000000 ## 3 5 1.0000000 0.0000000 ## 4 17 1.0000000 0.0000000

• M-step: Update the mixing probability vector, π1 = 0.9 + 1 + 1 + 0.9 + 1 + 0.692 6 = 0.915 π2 = 1 -π1 = 0.085, n <-nrow(Y) YP <-dplyr::left_join(Y, P_ZY, by = "target") YP

## # A tibble: 6 x 4 ## target spillover P_Y1 P_Y2 ## <dbl> <dbl> <dbl> <dbl> ## 1 3 2 0.9 0.1 ## 2 5 3 1 0 ## 3 17 2 1 0 ## 4 3 NA 0.9 0.1 ## 5 17 NA 1 0 ## 6 2 NA 0.692 0.308 pi <-c(sum(YP$P_Y1)/n, sum(YP$P_Y2)/n) pi ## [1] 0.91538462 0.08461538
and re-estimate the distribution for the target marker using the posterior probabilities as weights, keep the non-target marker at its initial value, The generative model for beads is an independent copy of Y | Z = 2.

P (Y i = 2 | Z i = 1) = 0.126 P (Y i = 2 | Z i = 2) = 2/3 P (Y i = 3 | Z i = 1) = 0.328 P (Y i = 3 | Z i = 2) = 1/3 P (Y i = 5 | Z i = 1) = 0.182 P (Y i = 17 | Z i = 1) = 0.364 F1

CFigure 2 :

 2 Figure 2: Three experiments testing our assumptions and sensitivy to bimodal bead distribution. For each experiment the top row are mean values over the entire range of the experimental setups, and the bottom row are density plots for three parameter settings to illustrate the generated distributions. Y is the distribution with spillover. Y | Z = 1 is the distribution without spillover. Y | Z = 2 is the spillover. mean(Y ) is the average of the distribution with spillover. mean(Y | Z = 1) is the average count without spillover. spillR mean(Y ) is the average count after correcting Y .

  <-spatstat.geom::ewcdf(Y$target, weights = YP$P_Y1) par(mfrow = c(1, 2)) plot(F1, main = "New CDF of Target", xlab = "y", ylab = "F1(y)") plot(F2, main = "CDF of Spillover", xlab = "y", ylab = "F2(y)") spillover probability estimate,P (spillover | Y i = 2) = 1 -0.915 • 0.126 0.915 • 0.126 + 0.085 • 2/3 = 0.33 P (spillover | Y i = 3) = 1 -0.915 • 0.328 0.915 • 0.328 + 0.085 • 1/3 = 0.09 P (spillover | Y i = 5) = 0 P (spillover | Y i = 17) = 0. P_ZY <-dplyr::mutate(P_YZ, P_Y1 = pi[1] * P_Y1, P_Y2 = pi[2] * P_Y2) P_ZY <-P_ZY /rowSums(P_ZY) P_ZY |> dplyr::mutate(target = c(2, 3, 5, 17)) |> dplyr::mutate(p_spillover = 1 -P_Y1) |> dplyr::select(target, p_spillover) Z = 2) ∼ Poisson(70 + τ ) (spillover component with shift)Y = (1 -I) • (Y | Z = 1) + I • (Y | Z = 2) (mixture).The generative model for beads is an independent copy of the unshifted Y | Z = 2 at τ = 0.Model MisspecificationGenerative model for real cells Y of this experiment:Z = 1) = (1 -M ) • T + M • S (target mixture component) (Y | Z = 2) = (1 -M ) S + M • T (spillover mixture component) Y = (1 -I) • (Y | Z = 1) + I • (Y | Z = 2) (mixture)The generative model for beads is an independent copy of Y | Z = 2.Bimodal SpilloverGenerative model for real cells Y of this experiment: Z = 2) = (1 -H) • (S | H = 0) + H • (S | H = 1) (spillover component) Y = (1 -I) • (Y | Z = 1) + I • (Y | Z = 2) (mixture)

Table 1 :

 1 Additional summaries of the data underlying Figure1C.

	correction	zeros or NA's mean
	none	31	4.45
	spillR	511	4.71
	other (CATALYST)	2162	2.94

https://github.com/marcoguazzini/spillR

https://github.com/ChristofSeiler/spillR_paper

Acknowledgments

We thank EuroBioC2022 for awarding Marco Guazzini a travel award to present a preliminary version of spillR in Heidelberg. We thank Antoine Chambaz for his feedback on an earlier draft that substantially improved the paper. Alexander G. Reisach received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 945332 .

<dbl> <dbl> ## 1 0.167 0.667 ## 2 0.333 0.333 ## 3 0.167 0 ## 4 0.333 0

We initialize the mixture probabilities with the discrete uniform, π1 = 0.9 π2 = 0.1.

pi <-c(0.9, 0.1)

Now, we update these initial values using the E and M-steps.

• E-step: Calculate the posterior probability for the true marker, and the spillover marker,


Saulnier Masson

Karim Daoues

John Measey

Anthony Herrel

The evolution of bite force and head morphology in scincid lizards: diet and habitat use as possible drivers

Keywords:

Diet and habitat use impose mechanical constraints that may impact head morphology and bite force. Skinks (Scincidae) comprise the largest family of lizards with ~1700 species currently described. They also show an important morphological and ecological diversity. Using phylogenetically informed analyses, we studied the interrelationships between ecology (diet, habitat use), head morphology, and bite force in these lizards. Our results show a strong link between body size, bite force, and diet, with herbivorous species being larger and biting harder than species from other dietary groups. Despite a lack of differences in body size and head morphology, omnivorous species bite harder than insectivorous species, likely to process the fibrous plant material that is part of their diet. Overall, lineages that evolved greater bite forces also showed an increased in relative head height allowing for more vertically oriented jaw muscles. Moreover, we find evidence for correlated evolution between bite force and head length: skinks that bite harder tend to have shorter jaws that likely provide a greater mechanical advantage when biting at the tip of the jaw. Surprisingly, habitat use does not appear to be correlated with morphological traits or bite force, but this needs to explored further.

Introduction

The vertebrate skull is a complex integrated system composed of numerous structural units (Herrel et al., 2007a. These units are forged by different biotic and abiotic interactions and are critical to several functions including defense (e.g., [START_REF] Cooper We | Incompletely protective refuges: selection and associated defences by a lizard, Cordylus cordylus (Squamata: Cordylidae)[END_REF], locomotion (e.g., [START_REF] Teodecki | Head dimorphism and burrowing speed in the African caecilian Schistometopum thomense (Amphibia: Gymnophiona)[END_REF][START_REF] Gans | Tetrapod limblessness: evolution and functional corollaries[END_REF], male-male combat (e.g., [START_REF] Huyghe | Morphology, performance and fighting capacity in male lizards, Gallotia galloti[END_REF][START_REF] Lappin | Gaping displays reveal and amplify a mechanically based index of weapon performance[END_REF], mating (e.g., Herrel et al., 1999b), drinking (e.g., [START_REF] Bels | Biomechanics of feeding in vertebrates[END_REF][START_REF] Cundall | Drinking in snakes: kinematic cycling and water transport[END_REF], and feeding (e.g., [START_REF] Wainwright | Predicting patterns of prey use from morphology of fishes[END_REF][START_REF] Cundall | Feeding in snakes. Feeding: form, function, and evolution in tetrapod vertebrates[END_REF]Schwenk, 2000). As such the skull plays a key role in many different ecological and behavioral contexts [START_REF] Hanken | The skull, Volume 3: Functional and evolutionary mechanisms[END_REF]. Moreover, the skull houses and protects the brain and the main sensory organs (e.g., eyes, vomeronasal apparatus). As these organs take up space within this integrated system there may be competing demands for space [START_REF] Barel | Towards a constructional morphology of cichlid fishes (Teleostei, Perciformes)[END_REF] and function (Herrel et al., 2001c) between these different structural units. Thus, the skull is the result of structural integration and functional compromises [START_REF] Wainwright | Predicting patterns of prey use from morphology of fishes[END_REF][START_REF] Vanhooydonck | Push and bite: trade-offs between burrowing and biting in a burrowing skink (Acontias percivali)[END_REF][START_REF] Corbin | Linkage and trade-off in trophic morphology and behavioural performance of birds[END_REF][START_REF] Edwards | Diving in head first: trade-offs between phenotypic traits and sand-diving predator escape strategy in Meroles desert lizards[END_REF][START_REF] Maestri | Diet, bite force and skull morphology in the generalist rodent morphotype[END_REF][START_REF] Watanabe | Ecomorphological diversification in squamates from conserved pattern of cranial integration[END_REF], which renders our understanding of the factors driving variation in cranial morphology complex.

The study of animal performance is important in this context as it allows variation in morphology to be linked to variation in ecology (Wainwright and Reilly, 1994;[START_REF] Irschick | Animal athletes: an ecological and evolutionary approach[END_REF]. Indeed, selection acts at the level of whole-organism performance rather than on the underlying structural components [START_REF] Arnold | Morphology, performance and fitness[END_REF][START_REF] Careau | Performance, personality, and energetics: correlation, causation, and mechanism[END_REF][START_REF] Lailvaux | The life history of whole-organism performance[END_REF]. Consequently, performance traits reflect both the behavior and the ecological context of a species [START_REF] Irschick | Intraspecific correlations among morphology, performance and habitat use within a green anole lizard (Anolis carolinensis) population[END_REF]Herrel et al., 2007a;[START_REF] Huyghe | Morphology, performance, behavior and ecology of three color morphs in males of the lizard Podarcis melisellensis[END_REF]. One of the most relevant performance traits associated with the cranial system is bite force [START_REF] Wainwright | Predicting patterns of prey use from morphology of fishes[END_REF][START_REF] Schluter | The ecology of adaptive radiation[END_REF][START_REF] Anderson | Bite force in vertebrates: opportunities and caveats for use of a nonpareil whole-animal performance measure[END_REF]. This performance trait is important in the acquisition of resources, in fighting, and in defense [START_REF] Herrel | Ecomorphology of the lizard feeding apparatus: a modelling approach[END_REF][START_REF] Erickson | The ontogeny of bite-force performance in American alligator (Alligator mississippiensis)[END_REF]. The use of novel resources and the conquest of new ecological niches by an animal are consequently often associated with morphological adaptations of the jaw system resulting in variation in bite force across species (Herrel et al., 2004b;[START_REF] Cattau | Rapid morphological change of a top predator with the invasion of a novel prey[END_REF]. Moreover, for some species, bite force has been demonstrated to be heritable [START_REF] Zablocki-Thomas | Maternal and genetic correlations between morphology and physical performance traits in a small captive primate, Microcebus murinus[END_REF] and to impact survival in the field [START_REF] Herrel | Do adult phenotypes reflect selection on juvenile performance? A comparative study on performance and morphology in lizards[END_REF], making it a key performance trait that is likely under selection.

Previous studies on lizards and other vertebrates have shown that the evolution towards higher bite force often goes hand in hand with the evolution of larger body size [START_REF] Aguirre | Ecomorphological analysis of trophic niche partitioning in a tropical savannah bat community[END_REF]Herrel et al., 2004a[START_REF] Herrel | Sexual dimorphism, body size, bite force and male mating success in tuatara[END_REF][START_REF] Chazeau | Proximate determinants of bite force capacity in the mouse lemur[END_REF]. However, irrespective of variation in body size, the evolution of a bigger head also promotes stronger bites due to the larger absolute jaw muscle volume that can be housed (Herrel et al., 2001a(Herrel et al., , 2001b[START_REF] Herrel | Ecological consequences of ontogenetic changes in head shape and bite performance in the Jamaican lizard Anolis lineatopus[END_REF]. Bite force has been suggested to determine diet in lizards with an increase in bite force often leading to a greater trophic diversity [START_REF] Wittorski | Proximate determinants of bite force in Anolis lizards[END_REF][START_REF] Taverne | From micro to macroevolution: drivers of shape variation in an island radiation of Podarcis lizards in the Adriatic[END_REF]. Durophagy and herbivory, for example, are dietary specialisations that require high bite force due to the mechanical resistance of these hard and fibrous food items (Herrel et al., 1999a;2004b;[START_REF] Schaerlaeken | Built to Bite: Feeding kinematics, bite forces, and head shape of a specialized durophagous lizard, Dracaena guianensis (Teiidae)[END_REF]. It has consequently been suggested that lizards eating such food items should have higher and wider heads with more massive cranial muscles than lizards that eat softer foods [START_REF] Herrel | Cranial morphology and bite force in Chamaeleolis lizards-adaptations to molluscivory?[END_REF][START_REF] Schaerlaeken | Built to Bite: Feeding kinematics, bite forces, and head shape of a specialized durophagous lizard, Dracaena guianensis (Teiidae)[END_REF].

On the other hand, habitat use can also generate mechanical limitations affecting head size and shape and thus indirectly also drive bite force evolution. In fossorial lizards, for example, the time needed to burrow, and the energetic cost of burrowing correlate with bite force given that both are determined by head width (e.g., [START_REF] Navas | Morphological and physiological specialization for digging in amphisbaenians, an ancient lineage of fossorial vertebrates[END_REF][START_REF] Barros | Head shape evolution in Gymnophthalmidae: does habitat use constrain the evolution of cranial design in fossorial lizards?[END_REF][START_REF] Vanhooydonck | Push and bite: trade-offs between burrowing and biting in a burrowing skink (Acontias percivali)[END_REF][START_REF] Guilloux | Tradeoffs between burrowing and biting force in fossorial scincid lizards?[END_REF]. In addition, crevice dwellers or climbing lizards may show reduced head and body heights to be able to exploit crevices and to avoid toppling backwards during climbing (Herrel et al., 2001b, c;[START_REF] Kohlsdorf | Head shape evolution in Tropidurinae lizards: does locomotion constrain diet?[END_REF]. Moreover, tree dwelling lizards typically have narrow heads allowing more stability when running on narrow branches (Herrel et al., 2001d;[START_REF] Kohlsdorf | Head shape evolution in Tropidurinae lizards: does locomotion constrain diet?[END_REF]. Consequently, head size and shape may be the result of trade-offs between feeding, habitat use, and locomotion and, as such, affect bite force [START_REF] Vanhooydonck | Push and bite: trade-offs between burrowing and biting in a burrowing skink (Acontias percivali)[END_REF]. Therefore, not only head size but also head shape is likely an important determinant of bite force (Herrel et al., 2007b;Fabre et al., 2014a, b).

In this study, we analyze the proximate determinants of bite force across 56 species of the family Scincidae. This group of lizards was chosen as it is the most species-rich lizard family characterized by an exceptional morphological and ecological diversity [START_REF] Chapple | Conservation status of the world's skinks (Scincidae): Taxonomic and geographic patterns in extinction risk[END_REF][START_REF] Uetz | The Reptile Database[END_REF]. Skinks show a great diversity in diet ranging from insectivorous to herbivorous and durophagous species.

Moreover, they can be found on all continents except Antarctica and exploit a wide variety of habitats ranging from arboreal over fossorial to terrestrial and even semi-aquatic. In addition to this tremendous ecological diversity, scincid lizards are also morphologically diverse, yet remain relatively poorly studied [START_REF] Greer | The genetic relationships of the scincid lizard genus Leiolopisma and its relatives[END_REF]William and Peterson, 1982;[START_REF] Paluh | Comparative skull anatomy of terrestrial and crevice-dwelling Trachylepis skinks (Squamata: Scincidae) with a survey of resources in scincid cranial osteology[END_REF][START_REF] Foster | Skink ecomorphology: forelimb and hind limb lengths, but not static stability, correlate with habitat use and demonstrate multiple solutions[END_REF].

Capitalizing on the morphological, ecological, and taxonomic diversity of skinks, we examine here the relationship between head morphology and bite force in the sample of skinks and assess whether this differs among lizards with distinct ecologies (diet, habitat use). We specifically predict that bite force will differ between diet groups with herbivorous species biting harder than insectivorous species as shown previously for other lizard groups (Herrel et al., 1999a;[START_REF] Vitt | History and the global ecology of squamate reptiles[END_REF]Herrel et al., 2004b;[START_REF] Metzger | Correlations between lizard cranial shape and diet: a quantitative, phylogenetically informed analysis[END_REF]. We further predict that fossorial species, climbers, and saxicolous species will differ in head size and shape and will have lower bite forces given the constraints on head size and shape in these habitats (Herrel et al., 2001d;[START_REF] Kohlsdorf | Head shape evolution in Tropidurinae lizards: does locomotion constrain diet?[END_REF][START_REF] Barros | Head shape evolution in Gymnophthalmidae: does habitat use constrain the evolution of cranial design in fossorial lizards?[END_REF][START_REF] Vanhooydonck | Push and bite: trade-offs between burrowing and biting in a burrowing skink (Acontias percivali)[END_REF][START_REF] Paluh | Comparative skull anatomy of terrestrial and crevice-dwelling Trachylepis skinks (Squamata: Scincidae) with a survey of resources in scincid cranial osteology[END_REF].

Material and Methods

Specimens

The sample consisted of data on head dimensions, body size, and bite force for 331 individuals across 56 species of skinks. Seven of these belong to the Acontinae, 16 to the Scincinae and 33 to the Lygosominae. We use binomial nomenclature as provided by the Reptile Database [START_REF] Uetz | The Reptile Database[END_REF].

The number of individuals sampled per species (Table 1) varies according to the availability of specimens for morphological and in vivo measurements.

Morphometrics

Five morphological measurements were taken for all specimens (Fig. 1, Table 1). We used the same measurements as described in [START_REF] Herrel | Cranial morphology and bite force in Chamaeleolis lizards-adaptations to molluscivory?[END_REF]. Snout-vent length (SVL) was measured from the tip of the snout to the posterior edge of the anal scale, head length (headl) from the back of the parietal bone to the tip of the upper jaw, head height (headh) was measured at the highest part of the head and posterior to the orbit, head width (headw) at the widest part of the head and at the level of jugal bone, and lower jaw length (lj) was measured from the back of the retroarticular process to the tip of lower jaw. The bones used in these measurements were easily identified by palpation. All measurements were taken using digital calipers (Mitutoyo CD-20DC, Kawasaki, Japan; precision: 0.01 mm), and were taken on the right side of the specimens.

Bite force

In vivo bite forces were taken in captivity or in the field. Only data for adults were used for this study and the number of individuals per species varied (Table 1). Bite force measurements were taken using an isometric Kistler force transducer (type 9203, range ±500 N; Kistler, Zurich, Switzerland) mounted on a purpose-built holder and connected to a Kistler charge amplifier (Type 5995 A, Kistler; see Herrel et al., 1999b). When the bite plates were placed between the jaws of the animals, prolonged and repeated biting typically resulted. The place of application of bite forces was standardized for all animals by metal stops that were mounted on the bite plates, thus assuring that animals always bit at the same position along the tooth row. Gape angle was standardized by moving the bite plates away from each other for larger animals resulting in animals biting at a gape angle of 30°. Surgical tape was wrapped around the bite plates to provide grip and to prevent damage to the teeth. Bite forces were measured five times for each animal and only the highest measurement was retained and considered to be the maximal bite force. The maximal raw bite force recorded from the amplifier was multiplied by 0.67 to correct for the lever arms of the set-up.

Diet

Species were classified into one of three dietary categories: herbivory (4 species), omnivory (7 species), and insectivory (45 species) (Table 2). Here, herbivorous species are considered to have a diet consisting of at least 70% fibrous plant material (i.e., all plant material except fruits and nectar), omnivorous species have a diet with 10% to 70% of plant material, and insectivorous species are those with less than 10% of plant material [START_REF] Cooper We | Distribution, extent, and evolution of plant consumption by lizards[END_REF]. As our study included a broad diversity of species, dietary data were not always available for each species. In that case the diet of a sister taxon of the same genus for which data were available was used. Sister taxa were used in only a few instances, and in all cases, diet did not show much variation within the genus. Moreover, we only extrapolated from other species in cases of insectivory to avoid inflating the number of independent origins of the less common diet types.

Habitat use

Skinks were classified into five habitat groups, in part based on [START_REF] Meiri | Traits of lizards of the world: Variation around a successful evolutionary design[END_REF]: arboreal (2 species), fossorial (22 species), saxicolous (3 species), semi-aquatic (3 species), and terrestrial (26 species) (Table 2). Arboreal species live in the vegetation and are rarely observed on the ground, unlike saxicolous species, which are as comfortable on the ground as in trees or rocks [START_REF] Ribeiro | Thermoregulatory behavior of the saxicolous lizard, Tropidurus torquatus (Squamata: Tropiduridae), in a rocky outcrop in Minas Gerais, Brazil[END_REF]. Semi-aquatic species spend significant amounts of time in water and are excellent swimmers, contrary to fossorial species which spend most of their time buried in the substrate. All other species were classified as terrestrial.

Statistical analysis

All statistical analyses were run in R (R Development Core Team, 2005). First, we calculated species means using only adults from the data set. Second, we pruned the time-calibrated phylogeny from Zheng and Wiens (2016) using the R packages ape [START_REF] Paradis | ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R[END_REF] and Geiger (Pennet et al., 2014) to include only the species in our data set (Fig. 2). To account for the non-independence of species due to shared ancestry [START_REF] Felsenstein | Phylogenies and the comparative method[END_REF], all statistical analyses were performed in a strict phylogenetic comparative framework.

To explore the effect of diet and habitat on overall body size (SVL), we performed a two-way phylogenetic analysis of variance (PhylANOVA) using a Brownian motion model followed by phylogenetic post-hoc tests using the aov.phylo function of geiger [START_REF] Pennell | geiger v2. 0: an expanded suite of methods for fitting macroevolutionary models to phylogenetic trees[END_REF] and the phylANOVA function of the phytools package [START_REF] Revell | Package 'phytools[END_REF].

To explore morphological differences according to diet and habitat, we performed multivariate phylogenetic analysis of covariance (phylogenetic MANOVA) using residual head dimensions. To do so all morphological variables were logarithmically transformed (log10) to fulfil assumptions of normality and homoscedasticity [START_REF] Kachigan | Multivariate statistical analysis: A conceptual introduction[END_REF][START_REF] Sokal | Biometry: The Principles and Practice of Statistics in Biological Research[END_REF]. Subsequently, the Log10-transformed head dimensions were regressed on Log10-transformed SVL using generalized least squares regressions with phylogenetic size correction (PGLS; [START_REF] Mao | Package 'pGLS': Generalized Least Square in comparative Phylogenetics[END_REF] and unstandardized residuals were extracted. These were then used as input for a phylogenetic MANOVA to test the effects of diet and habitat irrespective of variation in head size. Subsequently phylogenetic ANOVAs and post-hoc tests with phylolm function from phylolm package [START_REF] Ho | Package 'phylolm[END_REF] were run to better understand which variables drove the results and which groups differed from one another.

To determine which morphological traits explained variation in bite force, all morphological and bite force measures were log10-transformed and used as input for a phylogenetic stepwise multiple regression analysis with bite force as the dependent variable and the morphological traits as the independent variables. To do this, the phylostep function of phylolm package [START_REF] Ho | Package 'phylolm[END_REF] was used with Brownian Motion as the evolutionary model. This evolutionary model was selected by running the fitContinuous function from Geiger package (Pennet et al., 2014). The best regression model was selected based on the minimal AIC and variables in the model were selected using both forward and backward procedures. A phylogenetic regression was then using the phylolm function from the phylolm package [START_REF] Ho | Package 'phylolm[END_REF] to estimate the relationship between bite force and the morphological traits.

To investigate the effect of diet and habitat on bite force, we performed a two-way phylogenetic analysis of variance (PhylANOVA) with a Brownian motion model and followed by phylogenetic posthoc tests using the aov.phylo function of geiger [START_REF] Pennell | geiger v2. 0: an expanded suite of methods for fitting macroevolutionary models to phylogenetic trees[END_REF] and the phylANOVA function of the phytools package [START_REF] Revell | Package 'phytools[END_REF]. The post-hoc pairwise comparison tests taking phylogeny into account were performed using the Bonferroni method. As bite force is correlated with overall body size, we analyzed differences in bite force between diet and habitat groups using phylogenetic analysis of covariance (PhylANCOVA) with SVL as our co-variate. As no post-hoc pairwise comparisons are possible using this approach we extracted the unstandardized residuals from a bite force -SVL phylogenetic regression. These residuals (further referred to as residual bite force in the manuscript) were used as input for a PhylANOVA coupled to phylogenetic post-hoc pairwise comparison tests to test for differences between diet and habitat groups. Analyses were run separately for each diet and habitat use as our a priori predictions were different.

Results

The two-way PhylANOVA testing for differences in snout-vent length between species consuming different prey and utilizing different habitats showed significant differences between diet groups (F2, 49 = 20.14; P = 0.002), but not between species occupying different habitats (F4, 49 = 1.47; P = 0.58). The interaction between habitat and diet was also significant (F2, 47 = 7.11; P = 0.003). Phylogenetic posthoc tests indicated that herbivorous species were larger than omnivorous (P = 0.015) and insectivorous species (P = 0.003).

A phylogenetic MANOVA detected significant differences in head shape between diet groups (Wilks' Lambda = 0.36; F8, 100 = 8.23; P < 0.001). Subsequent univariate PhylANOVAs indicated that species were different for most variables with herbivorous species having relatively wider heads than insectivorous species. A second phylogenetic MANOVA further indicated significant differences in head shape among species from different habitats (Wilks' Lambda = 0.26; F16, 147.28 = 5.05; P = 0.018).

However, the univariate PhylANOVAs indicated no differences after Bonferroni correction. An inspection of the means suggests that fossorial skinks differ the most from skinks occupying other habitats and had the smallest morphological traits suggesting they have small heads for their body size.

A phylogenetic regression of bite force on snout-vent length was significant (P < 0.001; Fig 3) suggesting that larger species have higher bite forces. The phylogenetic stepwise multiple regression analysis with head dimensions as predictors of bite force retained a significant model explaining 78% of the variation of bite force (R 2 = 0.78; AIC(k=2) = 485). Head height (β = 9.68) was the best predictor and positively impacted bite force while head length (β = -3.02) negatively impacted bite force.

A two-way phylogenetic analysis of variance indicated absolute bite force differences between dietary groups (F2, 49 = 36.7; P < 0.001) but not habitat groups (F4, 49 = 2.20; P = 0.43). The interaction between diet and habitat was significant (F2, 47 = 9.37; P = 0.002). Bonferroni phylogenetic post-hoc tests indicated that herbivorous species were different from omnivorous and insectivorous species (all P < 0.01). Omnivorous species were also different from insectivorous species (P = 0.021). An inspection of the means showed that herbivorous species are characterized by higher bite forces followed by omnivorous, and insectivorous species (Figure 4).

A phylogenetic ANCOVA performed on the bite force indicated significant differences (F2, 52 = 11.23; P = 0.003) between diet, but not habitats groups (F4, 50 = 6.77; P = 0.078) when taking into account differences in snout-vent length. The phylogenetic ANOVA using residual bite forces indicated significant differences between diet groups (F2, 49 = 6.78; P = 0.499). The same PhylANOVA on habitat groups indicated, however, no significant differences (F4, 49 = 5.32; P = 0.12). Phylogenetic post-hoc tests showed that herbivorous species were different from omnivorous species (P = 0.036) and insectivorous species (P = 0.003) in residual bite force. Insectivorous species were also different from omnivorous species (P = 0.039). An inspection of the means showed that herbivorous species have the highest and insectivorous species the lowest relative bite forces.

Discussion

As the lizard cranial system is implicated in many functions (e.g., [START_REF] Baeckens | Fossorial and durophagous: implications of molluscivory for head size and bite capacity in a burrowing worm lizard[END_REF] that may require conflicting head morphologies (e.g., [START_REF] Paluh | Comparative skull anatomy of terrestrial and crevice-dwelling Trachylepis skinks (Squamata: Scincidae) with a survey of resources in scincid cranial osteology[END_REF], trade-offs can occur that constrain its evolution [START_REF] Vanhooydonck | Push and bite: trade-offs between burrowing and biting in a burrowing skink (Acontias percivali)[END_REF]. Based on data on head morphology and bite force for 56 species of ecologically diverse skinks, we examined the interrelationships between head morphology, bite force, and ecology to assess the ecological correlates of cranial design and performance.

Phylogenetic comparative analyses showed that cranial morphology and performance were determined by dietary specialisation but not habitat use.

Our results show that bite force is dependent on size with larger animals showing greater bite forces as expected based on scaling relationships. Herbivorous species stand out as being larger than insectivorous and omnivorous species, thus providing them with a functional advantage and allowing them to generate larger absolute bite forces (Herrel et al., 2001b;[START_REF] Shea | Diet of two species of bluetongue skink, Tiliqua multifasciata and Tiliqua occipitalis (Squamata: Scincidae)[END_REF]2014). The evolution towards a larger body size in herbivorous species could also be facilitated by reduced selection on agility which is essential for lizards hunting active prey like insects or small vertebrates [START_REF] Van Damme | Evolution of herbivory in lacertid lizards: effects of insularity and body size[END_REF].

Moreover, the ubiquitous nature of plants in the environment and a presumably low foraging cost may provide additional advantages promoting large body size in herbivorous lizards [START_REF] Pough | Lizard energetics and diet[END_REF]. Large body size provides several other advantages including a longer intestinal tract, higher thermal inertia, and a lower surface to volume ratio, and may thus allow herbivorous lizards to meet their energetic demands [START_REF] Pough | Lizard energetics and diet[END_REF][START_REF] Van Damme | Evolution of herbivory in lacertid lizards: effects of insularity and body size[END_REF]. Interestingly, no size difference was observed between omnivorous and insectivorous species. Omnivorous species are considered opportunistically herbivorous with much of their diet being based on arthropods. This may constrain their body size if agility is an important constraint during prey capture (Schwenk, 2000;[START_REF] Herrel | Lizard ecology: The evolutionary consequences of foraging mode[END_REF][START_REF] Hoppe | Less need for differentiation? Intestinal length of reptiles as compared to mammals[END_REF].

Our data further highlight that herbivorous species have relatively wider heads than insectivorous species. A wider head positively affects the maximum allowable muscle volume, and thus the crosssectional area of the jaw adductors [START_REF] Bowman | Morphological differentiation and adaptation in the Galapagos finches[END_REF][START_REF] Herrel | Evolution of bite force in Darwin's finches: a key role for head width[END_REF]. It is tempting to interpret this difference in head morphology as the result of natural selection. Indeed, fibrous plant material requires greater bite forces to be reduced before swallowing and thus a wider head may provide a solution to these physical constraints (Herrel and De Vree, 1999;Herrel et al., 1999c;2004b). Although omnivorous species did not differ from insectivorous species in overall head shape, they did show a greater bite force. This difference might be explained by the inclusion of fibrous plant material into the diet, as highlighted for other taxa (Herrel et al., 1999b;2001a;[START_REF] Metzger | Correlations between lizard cranial shape and diet: a quantitative, phylogenetically informed analysis[END_REF]. Indeed, it can be expected that the most demanding food item is the one driving the upper limit of bite force and as such omnivorous species can be expected to evolve a high enough force to allow them to reduce fibrous and tough plant matter. However, the generalist morphology in omnivorous species may be driven by conflicting demands imposed by dietary diversity (Herrel et al., 2004b). Thus, an increase in bite force but not head size could allow a non-specialist to switch to an omnivorous diet by allowing the inclusion of harder or tougher food items. We might have expected insectivorous species to differ morphologically from omnivorous species with a smaller head and longer snout [START_REF] Toyama | Evidence of adaptive evolution in the cranial morphology of Tropidurid lizards from coastal Peru[END_REF]. This lack of morphological differences suggests that omnivorous species have retained the ancestral morphology of insectivorous species. In a second step, a specialization towards herbivory or a durophagous diet may be allowed by developing a wider, taller head [START_REF] Herrel | Cranial morphology and bite force in Chamaeleolis lizards-adaptations to molluscivory?[END_REF][START_REF] Schaerlaeken | Built to Bite: Feeding kinematics, bite forces, and head shape of a specialized durophagous lizard, Dracaena guianensis (Teiidae)[END_REF][START_REF] Meyers | The evolution of bite force in horned lizards: the influence of dietary specialization[END_REF].

Interestingly, and unlike what has been observed for other taxa (e.g., [START_REF] Barros | Head shape evolution in Gymnophthalmidae: does habitat use constrain the evolution of cranial design in fossorial lizards?[END_REF][START_REF] Openshaw | Head shape evolution in monitor lizards (Varanus): interactions between extreme size disparity, phylogeny and ecology[END_REF], no differences in overall body size were observed between species occupying different habitats. This suggests that habitat use may not impose strong constraints on size in skinks. Surprisingly, our results also indicated no differences in head shape and bite force between habitat groups other than a tendency for burrowers to have narrower heads. The tendency for burrowing species to differ in head morphology is not unexpected as a fossorial lifestyle may impose significant energetic constraints on head size [START_REF] Teodecki | Head dimorphism and burrowing speed in the African caecilian Schistometopum thomense (Amphibia: Gymnophiona)[END_REF][START_REF] Navas | Morphological and physiological specialization for digging in amphisbaenians, an ancient lineage of fossorial vertebrates[END_REF]. The lack of strong common morphological and functional specializations in species occupying different habitats is probably the result of functional compromises (e.g., fossorial and durophagous trade-offs; [START_REF] Baeckens | Fossorial and durophagous: implications of molluscivory for head size and bite capacity in a burrowing worm lizard[END_REF] or of independent evolutionary trajectories of different habitat specialists and needs to be investigated further. More, since habitat can impact the width or height of the body as well as limb morphology, further analyzes of body shape are needed to better understand the constraints imposed by different habitat types on morphology in skinks (e.g., [START_REF] Pounds | Ecomorphology, locomotion, and microhabitat structure: patterns in a tropical mainland Anolis community[END_REF][START_REF] Goodman | Life on the rocks: habitat use drives morphological and performance evolution in lizards[END_REF]Herrel et al., 2008;[START_REF] Meiri | Evolution and ecology of lizard body sizes[END_REF][START_REF] Pafilis | 6th Symposium on the Lacertids of the Mediterranean Basin[END_REF].

Head height and head length were the primary determinants of bite force in the species of skinks included in our study. A taller head, just like a wider head, likely increases the space available for jaw adductors which may drive the observed increase in bite force [START_REF] Bowman | Morphological differentiation and adaptation in the Galapagos finches[END_REF]Herrel et al., 2001b;[START_REF] Shea | Diet of two species of bluetongue skink, Tiliqua multifasciata and Tiliqua occipitalis (Squamata: Scincidae)[END_REF]. Moreover, a taller head may also allow for more vertically placed jaw adductors providing them with a greater moment arm and a more efficient conversion ratio of muscle force into bite force (Herrel et al., 2001b). Conversely, head length had a negative effect on bite force (see also [START_REF] Westneat | Evolution of levers and linkages in the feeding mechanisms of fishes[END_REF]).

An increase in head length involves an increase in the jaw outlever and given that animals were made to bite at the tips of the jaws this can be expected to negatively impact bite force (Herrel et al., 2001b;2007a). Indeed, the jaws of lizards are in the form of a third-order lever, with the lower jaw being the lever arm [START_REF] Kerr | Introductory biomechanics[END_REF][START_REF] Cox | The jaw is a second-class lever in Pedetes capensis (Rodentia: Pedetidae)[END_REF]. Thus, the variation in the length of the head will influence the distance between the input force and the output force of the lever [START_REF] Wainwright | Predicting patterns of prey use from morphology of fishes[END_REF]Fabre et al., 2014a). However, to better understand how head length may impact bite force a more indepth analysis of the jaw in-and outlevers in relation to bite force is needed.

Conclusion

Our data demonstrate a strong link between body size, bite force, and diet with herbivorous species being larger and biting harder in absolute terms. Moreover, in Scincidae the evolution towards more powerful bites goes hand in hand with an increase in head height but a decrease and head length.

Habitat use did not appear to be an important driver of head shape or bite force in the species included in our data set, yet this remains to be explored further using a broader and more comprehensive sample of species.

Zheng Y, Wiens JJ. 2016. Combining phylogenomic and supermatrix approaches, and a timecalibrated phylogeny for squamate reptiles (lizards and snakes) based on 52 genes and 4162 species. Molecular phylogenetics and evolution 94: 537-547. [START_REF] Fitzsimons | The lizards of South Africa: Family 4. Scincidae[END_REF][START_REF] Branch | Reptiles and amphibians of the Addo Elephant National Park[END_REF][START_REF] Bates | Atlas and red list of the reptiles of South Africa, Lesotho and Swaziland[END_REF] Scelotes sexlineatus Scincinae Insectivorous Fossorial [START_REF] Fitzsimons | The lizards of South Africa: Family 4. Scincidae[END_REF][START_REF] Branch | Reptiles and amphibians of the Addo Elephant National Park[END_REF][START_REF] Bates | Atlas and red list of the reptiles of South Africa, Lesotho and Swaziland[END_REF]

Figure legends

Figure 1 :

 1 Figure 1: Illustration of the measurements taken on the heads of the lizards, illustrated on a picture of the head of Eumeces schneideri.

Figure 2 :

 2 Figure 2: Time-calibrated phylogeny representing the relationships between the species included in this study modified from Zheng and Wiens (2016). The family, diet, habitat, and body form of each species are indicated.

Figure 3 :

 3 Figure 3: Scatter plot of species mean bite force against species mean snout-vent length. The diet and the habitat of species are indicated. The species depicted on the graph are: Corucia zebrata, Tiliqua occipitalis, Acontias meleagris, Typhlosaurus vermis, Scincus mitranus, Trachylepis quinquetaeniata and Ablepharus kitaibelii. The shaded area represents 95% confidence intervals.Note that x-and y-axes are on a logarithmic scale.

Figure 4 :

 4 Figure 4: Plot of bite force against snout-vent length. The regression lines for each diet group are represented. Red dots: insectivorous species; blue dots: omnivorous species; green dots: herbivorous species. The shaded areas represent 95% confidence intervals. Note the log-axes.

Table 1 :

 1 Morphological traits and bite force of the species included in our study.

	586
				Head length	Head width	Head height	lower jawl	Bite
	Species	n	svl (mm)	(mm)	(mm)	(mm)	length (mm)	force (N)
								0.3 ±
	Ablepharus kitaibelii	13	40.94 ± 2.43	6.44 ± 0.45	3.62 ± 0.28	2.55 ± 0.22	5.85 ± 0.57	0.14
								1.16 ±
	Acontias kgalagadi	5	116.12 ± 12.1	7.27 ± 0.49	3.6 ± 0.2	3.23 ± 0.86	5.7 ± 0.34	0.4
								0.77 ±
	Acontias litoralis	42 116.62 ± 11.72 6.22 ± 0.86	3.13 ± 0.24	2.66 ± 0.18	5.41 ± 0.62	0.19
	Acontias meleagris	6 194.24 ± 16.34 10.43 ± 1.08	5.81 ± 0.9	4.56 ± 0.55	7.67 ± 1.16	9 ± 1.77
								10.08 ±
	Acontias percivali	13 224.3 ± 10.04 14.26 ± 0.15	8 ± 0.06	6.98 ± 0.08	13.25 ± 0.44	2.62
	Bellatorias frerei	1	188.36	31.56	22.48	17.9	35.7	81.47
								8.51 ±
	Chalcides ocellatus	13 96.11 ± 24.98 13.46 ± 2.84	9.03 ± 2.24	7.7 ± 2	14.71 ± 3.18	4.19
								4.12 ±
	Chalcides sepsoides	17	83.41 ± 5.73	9.9 ± 0.48	5.93 ± 0.33	5.05 ± 0.24	9.97 ± 0.71	0.75
	Chalcides sphenopsiformis	1	83	7.99	4.97	4.25	7.51	3.16
								7.49 ±
	Cophoscincopus greeri	4	65.73 ± 2.09	12.93 ± 0.55	8.19 ± 0.66	6.98 ± 0.56	12.14 ± 1	0.79
								206.85 ±
	Corucia zebrata	2	275 ± 7.07	50.6 ± 4.53	48.88 ± 2.37 36.03 ± 3.01	57.1 ± 1.63	94.36
								7.86 ±
	Ctenotus uber	3	86.26 ± 23.66 15.08 ± 3.61	9.99 ± 2.41	8.81 ± 2.54	17.5 ± 4.19	4.19
	Cyclodomorphus gerrardii	1	70.92	14.9	10.15	7.34	17.15	16.683
								9.74 ±
	Cyclodomorphus michaeli	2 122.69 ± 14.11 15.49 ± 2.12	9.41 ± 1.07	9.12 ± 0.04	19.49 ± 2.7	1.15
								16.66 ±
	Eulamprus heatwolei	5	109.17 ± 23.2	21.37 ± 3.6	13.7 ± 2.9	11.24 ± 2.26	25.96 ± 5.23	7.24
								53.27 ±
	Eumeces schneiderii	3	122.01 ± 22.8 24.61 ± 1.57 16.16 ± 2.98 13.97 ± 2.53	27.48 ± 3.02	6.25
	Isopachys gyldenstolpei	1	180	10.54	6.6	5.82	8.77	5.89
								2.36 ±
	Leptosiaphos kilimensis	5	69.81 ± 14.29 10.26 ± 1.14	5.63 ± 0.84	4.18 ± 0.61	10.5 ± 1.55	1.87
								2.51 ±
	Lerista edwardsae	16	81.03 ± 7.34	7.99 ± 0.45	4.91 ± 0.35	4.18 ± 0.34	7.4 ± 0.87	0.63
								41.07 ±
	Mochlus fernandi	2 124.75 ± 11.92 22.09 ± 1.94 15.96 ± 2.21 12.19 ±2.04	25.36 ± 3.46	0.95
								8.84 ±
	Mochlus sundevallii	15 87.31 ± 29.36 12.73 ± 2.84	7.87 ± 1.87	6.37 ± 1.73	12.95 ± 3.08	5.06
	Morethia butleri	1	45.92	8.34	4.99	3.98	8.92	1.07
	Panaspis togoensis	1	35.59	6.89	3.4	2.2	6.95	0.33
	Panaspis wahlbergii	1	35.59	6.89	3.4	2.2	6.95	0.33
								6.58 ±
	Plestiodon fasciatus	4	67.72 ± 1.68	13.8 ± 1.32	9.31 ± 0.95	6.33 ± 0.76	15.72 ± 1.04	1.51
	Plestiodon gilberti	1	89.92	16.77	11.32	8.39	18.86	10.19
								17.86 ±
	Plestiodon laticeps	17 99.87 ± 11.14 19.59 ± 2.56 15.09 ± 2.84	9.75 ± 1.77	22.43 ± 3.01	4.26
								1.61 ±
	Plestiodon tetragrammus	2	58.72 ± 3.23	9.42 ± 0.25	6.38 ± 0.06	5.33 ± 0.31	11.36 ± 0.99	0.33
	Pygomeles braconnieri	1	165	13.31	8.17	6.65	12.33	9.36
								1.38 ±
	Scelotes bipes	9	75.47 ± 11.8	6.2 ± 0.62	3.18 ± 0.61	2.7 ± 0.32	5.53 ± 0.62	0.71
	Scelotes limpopoensis	1	61.07	7.3	3.81	2.81	7.32	0.51

  Table entries are means ± standard deviations. Morphological measurements are in mm, bite forces

	587
	in Newton. n: number of specimens measured.

Table 2 :

 2 Ecological traits for each species included in the study and references.

	589
	Species	Family	Diet	Habitat	References
	Ablepharus kitaibelii	Lygosominae Insectivorous Terrestrial Meiri, 2018
	Acontias kgalagadi	Acontinae	Insectivorous Fossorial	Meiri, 2018
	Acontias litoralis	Acontinae	Insectivorous Fossorial	Meiri, 2018
	Acontias meleagris	Acontinae	Insectivorous Fossorial	Meiri, 2018
	Acontias percivali	Acontinae	Insectivorous Fossorial	Meiri, 2018
	Bellatorias frerei	Lygosominae Omnivorous	Terrestrial Meiri, 2018
					Andrews et al., 1987; Carretero et
	Chalcides ocellatus	Scincinae	Insectivorous Fossorial	al., 2010
	Chalcides sepsoides	Scincinae	Omnivorous	Fossorial	Meiri, 2018
					Andrews and Kenney, 1990; Attum
	Chalcides sphenopsiformis Scincinae	Insectivorous Fossorial	et al., 2004; Meiri, 2018
				Semi-
	Cophoscincopus greeri	Lygosominae Insectivorous	aquatic	Meiri, 2018
	Corucia zebrata	Lygosominae Herbivorous	Arboreal	Meiri, 2018
	Ctenotus uber	Lygosominae Omnivorous	Terrestrial Meiri, 2018
	Cyclodomorphus gerrardii Lygosominae Insectivorous Terrestrial Meiri, 2018
	Cyclodomorphus michaeli Lygosominae Insectivorous Terrestrial Meiri, 2018
				Semi-
	Eulamprus heatwolei	Lygosominae Insectivorous	aquatic	Meiri, 2018
	Eumeces schneiderii	Scincinae	Omnivorous	Terrestrial Meiri, 2018
					Das, 2010; Chan-ard et al., 2015;
	Isopachys gyldenstolpei	Lygosominae Insectivorous Fossorial	Camaiti et al., 2022
	Leptosiaphos kilimensis	Lygosominae Insectivorous Terrestrial Meiri, 2018
	Lerista edwardsae	Lygosominae Insectivorous Fossorial	Meiri, 2018
	Mochlus fernandi	Lygosominae Insectivorous Terrestrial Meiri, 2018
	Mochlus sundevallii	Lygosominae Insectivorous Terrestrial Meiri, 2018
	Morethia butleri	Lygosominae Insectivorous Terrestrial Meiri, 2018
	Panaspis togoensis	Lygosominae Insectivorous Terrestrial Meiri, 2018
					Razzetti and Msuya, 2002; Ceriaco et
	Panaspis wahlbergii	Lygosominae Insectivorous Terrestrial	al., 2018
	Plestiodon fasciatus	Scincinae	Insectivorous Terrestrial Meiri, 2018
	Plestiodon gilberti	Scincinae	Insectivorous Terrestrial Meiri, 2018
	Plestiodon laticeps	Scincinae	Insectivorous Terrestrial Meiri, 2018
	Plestiodon tetragrammus Scincinae	Insectivorous Terrestrial Meiri, 2018
	Pygomeles braconnieri	Scincinae	Insectivorous Fossorial	Glaw and Vences, 2007; Meiri, 2018
	Scelotes bipes	Scincinae	Insectivorous Fossorial	Meiri, 2018
	Scelotes limpopoensis	Scincinae	Insectivorous Fossorial	Meiri, 2018
	Scelotes montispectus	Scincinae	Insectivorous Fossorial

  2014

	Scincus scincus	Scincinae	Omnivorous	Fossorial	Meiri,
	Tiliqua multifasciata	Lygosominae Herbivorous	Terrestrial Meiri,
	Tiliqua occipitalis	Lygosominae Herbivorous	Terrestrial Shea, 2006; Swan et al., 2017
	Tiliqua rugosa	Lygosominae Herbivorous	Terrestrial Meiri,
	Tiliqua scincoides	Lygosominae Omnivorous	Terrestrial Meiri,
	Trachylepis capensis	Lygosominae Insectivorous Terrestrial Meiri,
	Trachylepis
	homalocephala	Lygosominae Insectivorous Saxicolous Meiri,
	Trachylepis punctatissima Lygosominae Insectivorous Terrestrial Meiri,
	Trachylepis				Spawls et al., 2001; Meiri, 2018;
	quinquetaeniata	Lygosominae Insectivorous Saxicolous	Dendi et al., 2019
	Trachylepis spilogaster	Lygosominae Insectivorous Arboreal	Meiri,
	Trachylepis striata	Lygosominae Insectivorous Terrestrial Meiri,
	Trachylepis sulcata	Lygosominae Insectivorous Saxicolous Meiri,
	Trachylepis varia	Lygosominae Insectivorous Terrestrial Meiri,
	Trachylepis variegata	Lygosominae Insectivorous Terrestrial Meiri,
	Tribolonotus gracilis	Lygosominae Insectivorous Terrestrial Meiri,
	Tribolonotus
	novaeguineae	Lygosominae Insectivorous Terrestrial Meiri.
				Semi-	Hikida et al., 2003; Chuaynkern et
	Tropidophorus baconi	Lygosominae Insectivorous	aquatic	al., 2014
	Typhlosaurus caecus	Acontinae	Insectivorous Fossorial	Meiri,
	Typhlosaurus lomiae	Acontinae	Insectivorous Fossorial	Meiri,
	Typhlosaurus vermis	Acontinae	Insectivorous Fossorial	Meiri,
	590
	591
	Scincella cherriei	Lygosominae Insectivorous Terrestrial Meiri, 2018
	Scincopus fasciatus	Scincinae	Insectivorous Fossorial	Meiri, 2018
	Scincus mitranus	Scincinae	Omnivorous	Fossorial	Meiri, 2018

Table 4 :

 4 Bonferroni post-hoc results testing for differences between for diet categories.Significant results are in bold. SE; Standard Errors

		Herbivorous	Insectivorous	Omnivorous
	Snout-vent length
	Herbivorous		0.003	0.015
	Insectivorous	0.003		0.534
	Omnivorous	0.015	0.534
	Bite force
	Herbivorous		0.003	0.003
	Insectivorous	0.003		0.021
	Omnivorous	0.003	0.021
	Revaltive bite force
	Herbivorous		0.594	1
	Insectivorous	0.594		0.15
	Omnivorous	1	0.15
	Head length
	Herbivorous		0.222	1
	Insectivorous	0.222		0.282
	Omnivorous	1	0.282
	Head height
	Herbivorous		0.075	1
	Insectivorous	0.075		0.126
	Omnivorous	1	0.126
	Head width
	Herbivorous		0.048	0.528
	Insectivorous	0.048		0.396
	Omnivorous	0.528	0.396
	Lower jaw length
	Herbivorous		0.33	1
	Insectivorous	0.33		0.477
	Omnivorous	1	0.477
	Significant results are in bold.

Acknowledgements

We would like to thank Simon Baeckens and an anonymous reviewer for helpful and constructive comments on an early version of the manuscript; Krystal Tolley for help in the field and for collecting animals used in the study, the South African National Research Foundation (Key International Science Capacity Fund Program), the Partenariat Hubert Curien -National Research Foundation Protea, the South African National Biodiversity Institute, the City of Cape Town, and Oppenheimer and Sons for financial and logistical support. Field work in Australia was conducted under permit number E24474-1.

Data availability

All data are available in table one of the manuscript.


Nicolas Anciaux
email: nicolas.anciaux@inria.fr

Sabine Frittella
email: sabine.frittella@insa-cvl.fr

Baptiste Joffroy
email: baptiste.joffroy@insa-cvl.fr

Benjamin Nguyen
email: benjamin.nguyen@insa-cvl.fr

Demo : Data Minimization and Informed Consent in Administrative Forms

Keywords: Security and privacy → Privacy protections Privacy, Data minimization, GDPR, Informed consent

This article proposes a demonstration implementing the data minimization privacy principle, focusing on reducing data collected by government administrations through forms. Data minimization is defined in many privacy regulations worldwide, but has not seen extensive real-world application. We propose a model based on logic and game theory and show that it is possible to create a practical and efficient solution for a real French welfare benefit case.

INTRODUCTION

Overview. In this demonstration, we propose a tool to assist users with the problem of the minimization of personal data collected in the context of administrative forms. Data minimization is a fundamental privacy concept defined in many privacy regulations (such as the European General Data Protection Regulation or GDPR [START_REF] Council | Regulation EU 2016/679 of the European Parliament and of the Council[END_REF]) that states that only the data necessary to take a decision must be collected and stored. Data minimization has yet to see any kind of large scale (non trivial) implementation. We have recently proposed in [START_REF] Anciaux | A new PET for data collection via forms with data minimization, full accuracy and informed consent[END_REF] a model, based on propositional logic and game theory, to capture the general problem of data minimization. We propose to demonstrate the applicability of this approach on real forms.

General Context. Forms are filled by individuals who wish to apply for welfare benefits (e.g. in France). The volume of personal data thus collected is massive and concerns millions of individuals. For example, hundreds of different types of aids are offered in France at local or national level, to meet different social needs or provide incentives, with a very large number of beneficiaries. This includes aid for family, education, savings, access to housing, health, transport, mobility, energy saving, legal aid, etc. Governmental portals offer simulators so that citizens can identify the aids they can apply for. For example, aide-sociale.fr offers a rather holistic simulator (more than 1000 aids); the governmental portal mesdroitssociaux.gouv.fr is targeted on the most solicited aids (58 welfare aids) and also offers a simulator for young adults (887 aids). In order to be obtained, each proposed aid requires the beneficiary to fill out a form to collect their personal data showing that they meet the eligibility criteria and calibrating the aid. The allocation decision is made based on the collected data, which leads to the processing of millions of forms per year, with several weeks of instruction by the administration. For example, the public data set provided by the French government indicates that 6 million forms per year were processed on average over the period 2016-2020 for only 10 family welfare benefits. As another example, the social aid benefit related to "complementary health coverage" (which is one of those we propose to demonstrate) concerns 7.19 million beneficiaries in 2022 (see annual report of the complementary health insurance on p.10), all of whom had to send in their corresponding form to activate and subsequently renew the aid each year (see same report at bottom of p.8).

Objective. Given these staggering numbers, an effective solution to minimize data collected via forms, while properly informing the individual and obtaining their informed consent to the collection process, may lead to significant resource gains and more rigorous compliance with legal principles [START_REF] Determann | The california privacy rights act of 2020: A broad and complex data processing regulation that applies to businesses worldwide[END_REF][START_REF] Council | Regulation EU 2016/679 of the European Parliament and of the Council[END_REF]. We will show in this demonstration a PET (Privacy Enhancing Technology) which analyzes a given user's data, proposes all possible data minimizations, and quantifies the information that will subsequently be sent by the user, in order to let them make an informed decision on which minimization to choose.

DATA MINIMIZATION RELATED WORK

Data minimization prior to processing has long been considered an unsolvable problem [START_REF] Tene | Privacy in the age of big data: a time for big decisions[END_REF], with an overly negative influence on the quality of the processing [START_REF] Custers | Big data and data reuse: a taxonomy of data reuse for balancing big data benefits and personal data protection[END_REF], sometimes in contradiction with other legal obligations such as discrimination prevention [START_REF] Galdon Clavell | Auditing algorithms: On lessons learned and the risks of data minimization[END_REF], fairness [START_REF] Ramadan | A semi-automated bpmn-based framework for detecting conflicts between security, data-minimization, and fairness requirements[END_REF] or fraud detection [START_REF] Găbudeanu | Privacy intrusiveness in financial-banking fraud detection[END_REF]. Implementation is deemed too complicated for developers [START_REF] Alhazmi | I'm all ears! listening to software developers on putting gdpr principles into software development practice[END_REF][START_REF] Senarath | Understanding user privacy expectations: A software developer's perspective[END_REF], which has led to an overall limited adoption and compliance with this principle. Some barriers were removed by the introduction of a formalism clarifying the general objectives of a data minimizer [START_REF] Antignac | Data minimisation: a language-based approach[END_REF]. Indeed, recent contributions show that such a formalism can be applied in the context of using machine learning to predict recommendations/personalization through performance metrics [START_REF] Biega | Operationalizing the legal principle of data minimization for personalization[END_REF][START_REF] Shabanian | Learning to limit data collection via scaling laws: Data minimization compliance in practice[END_REF], while minimizing the data taken as input and still producing a useful model as output, even though its quality is slightly reduced. In [START_REF] Anciaux | A new PET for data collection via forms with data minimization, full accuracy and informed consent[END_REF], we study the case of administrative forms, where such a trade-off between data minimization (or privacy) and performance (utility of the processing) does not apply. Indeed the full range of benefits due to individuals must be offered by virtue of law: any metric that reduces the legitimate benefits due to a user in exchange for better data minimization would not be acceptable. Moreover, using logic and game theory, we show how to bring meaningful information to the individual in order to obtain their informed consent on the collection process, which is also required by regulations. This demonstration is based on this last result. We explain our approach on an example in the following section, which will be used in the demonstration.

A DATA MINIMIZATION EXAMPLE

The RSA benefit. In this section, we walk the reader through our approach, using the "active solidarity income" (in French: "revenu de solidarité active", abbreviated as RSA) welfare benefit form. The volume of RSA requests is around 2 million forms processed per year 1 for a total spending of around 12 billion euros, which is the highest welfare benefit (in term of euros spent) in France. Note that our demonstration will feature the RSA, but also the complementary health coverage aids which is the welfare benefit with the most forms processed annually in France (over 7 million forms) for a total spending of around 2.8 billion euros, as detailed in [START_REF] Anciaux | A new PET for data collection via forms with data minimization, full accuracy and informed consent[END_REF]. We believe these examples show that our approach can be used to manage real large scale use cases, and that there is a real impact on 1) the minimization of data sent and 2) the quantification of the data sent, thus improving the users information with a view to consent.

General approach. We consider an automatic decision making system, which processes user data, collected as a form of attribute/value pairs (e.g. 𝑎𝑔𝑒 = 25 or 𝑝𝑟𝑒𝑔𝑛𝑎𝑛𝑡 = 𝑓 𝑎𝑙𝑠𝑒). We assume that the automatic decision making system uses logical rules, which we write in disjunctive normal form (DNF) (e.g. a disjunction of conjunctions, each conjunction representing one possibility of obtaining a benefit). This is the case of many classical AI approaches such as decision trees, which benefit from explainability. By considering the values of all a given user's attributes, we deduce the predicates triggered by the user.

Building the RSA form. The RSA example was built using data from its official description 2 . 17 predicates are used, which are a direct translation of the official texts, negations included : 1 See the Court of Auditors website, and for more details its recent report: "The Active Solidarity Income", Evaluation of a public policy, Jan. 2022 (link to the report). 2 See here for the official description (in French). Eligibility criteria and benefits for RSA are outlined on the French "Family Allowance Fund". We directly translated the rules to a propositional formula, and then converted the formula to DNF using the Logical Expression Converter program (Haskell code available here). We then use a parser to transform this expression, and combine it with 𝑅 𝐴𝐷𝐷 to produce the DNF format used by our application (which is a DNF version of the DIMACS CNF format 3 ). The RSA.dnf file is available here. It is composed of the 17 predicates and a disjunction of 31 conjunctions of up to 12 different predicates. An example rule to obtain the benefit (e.g. the first of the 31 conjunctions) states that one needs to be French (𝑝 1 ), be over 25 years old (𝑝 5 ), live in France (𝑝 10 ), be a single parent (𝑝 11 ), and be neither conducting an internship (𝑝 14 ), a sabbatical (𝑝 16 ) or non-active (𝑝 17 ). The final rule is obtained by adding the predicates triggered in 𝑅 𝐴𝐷𝐷 , thus the final rule is : 𝑝 1 ∧𝑝 2 ∧¬𝑝 4 ∧𝑝 5 ∧¬𝑝 8 ∧𝑝 10 ∧𝑝 11 ∧𝑝 14 ∧𝑝 16 ∧𝑝 17 .

𝑝 1 : nationality = FRENCH
Computing the Minimal Accurate Subvaluations. We precompute (offline) the Minimal Accurate Subvaluations (MAS) for all possible forms that can obtain the benefit. MAS are a subset of the 17 predicates that correspond to enough information to grant the benefit to a user, even if the user had multiple ways of obtaining the benefit. These subvaluations are minimal in the sense that no predicate can be removed while still granting the benefit to the user. In this example, there are 1296 cases to consider, and 24 MAS. One example of a MAS would be noted 𝑚 1 = _1001__0_11___111 where a 0 or 1 in a given position means the predicate must be revealed, along with its value true (1) or false (0) and _ (blank) means the predicate is not sent.

Linking MAS to forms. We next build a bipartite graph linking each MAS to all the forms that it is a subvaluation of. It is possible that a MAS is linked to several forms, and it is possible that a form is linked to several MAS. If a form is only linked to one MAS, then Choosing the MAS using game theory. In order to choose which MAS to send, we need to define a cost function, for instance the number of forms sending the MAS. We note this payoff function 𝑃𝑂 𝑆𝑀 (𝑚) (Same MAS) where 𝑚 is a MAS. We use game theory to compute the correct value of 𝑃𝑂 𝑆𝑀 . The game is to decide which MAS each player (a form) will send to the administration, and the objective of the game is to send the MAS with the highest 𝑃𝑂 𝑆𝑀 value. This is a one round game, all players announcing simultaneously their move. We show in [START_REF] Anciaux | A new PET for data collection via forms with data minimization, full accuracy and informed consent[END_REF], Theorem 4.6 that this game admits a Nash equilibrium (thus the moves of rational players can be accurately predicted and computed), under a hypothesis on lexicographical ordering of the MAS : all players with only one choice play this choice, all players with equivalent choices will play the first MAS in lexicographical order. Note that each player reduces the value of 𝑃𝑂 𝑆𝑀 (𝑚) for all the 𝑚 that the player could have played, but in fact did not.

PET service for data collection using forms

Informed Data Minimization

Example. Consider a user with the form (as in the demo video) : 𝑓 = 1100110111000111 (a 0 or 1 in a position means a false or true value for the 𝑖 𝑡ℎ predicate, bear in mind some are negative such as 𝑝 15 = 𝑡𝑟𝑢𝑒, which means the person is not on leave). This user has the choice between 2 different MAS : 𝑚 1 = _1001__0_11___111 which corresponds to 32 other different forms or 𝑚 2 = _100_1___11___111 which corresponds to 128 other different forms. Note that the number of different forms is less or equal to 2 𝑏 where 𝑏 is the number of blank predicates. Our informed consent minimizer thus informs the user of this possible choice, and lets the user decide, suggesting the better choice of 𝑚 2

DEMONSTRATION OUTLINE

Our demonstration is divided into 3 parts. The overall architecture is presented in Fig. 1.

(1) Presentation/creation of the decision process. We have already analyzed 2 French administrative welfare benefit procedures (RSA and Welfare benefits). It is also possible for participants to create their own decision process rules. (2) Processing of the decision process rules (PET service for data collection). This demonstrates the feasibility of our approach on real world examples. Processing of the examples takes up to a few minutes on a regular laptop. This process constructs the MAS then generates the bipartite graph (Algo. 1) and computes the payoff function using game theory (Algo. 2). (3) Informed consent for users (Data collection service/GUI).

Participants are invited to fill in a form. The form is then processed using the output of the previous service. Users are suggested minimized forms, and are informed of the corresponding payoff functions (i.e. number of other players who share the same "blanked" form).

CONCLUSION

We introduced a novel Privacy Enhancing Technology for data minimization in administrative forms, informing eligible users about data removal and its impact. This confirms data minimization's practicality in decision rule contexts.

𝑝

Figure 1 :

 1 Figure 1: Informed Consent Data Minimization Architecture

  2 : nationality = EUROPEAN COUNTRY 𝑝 3 : nationality = SWISS 𝑝 4 : valid residence permit for at least 5 years = TRUE 𝑝 5 : age ≥ 25 𝑝 6 : pregnant = TRUE 𝑝 7 : children > 0 𝑝 8 : 18 ≤ age < 25 𝑝 9 : working = TRUE 𝑝 10 : living in France = TRUE 𝑝 11 : single parent = TRUE 𝑝 12 : High school student = FALSE 𝑝 13 : Higher education student = FALSE 𝑝 14 : Intern = FALSE 𝑝 15 : on leave = FALSE 𝑝 16 : sabbatical = FALSE 𝑝 17 : non active status = FALSE As some predicates are linked to each other, we add some consistency rules, noted 𝑅 𝐴𝐷𝐷 , as follow : 𝑝 1 → 𝑝 2 (if you are French you are from Europe) 𝑝 1 → ¬𝑝 4 (if you are French you do not have a green card) 𝑝 4 → ¬𝑝 1 (if you have a green card you are not french) 𝑝 5 → ¬𝑝 8 (if you are over 25, then you are not 18-25) 𝑝 8 → ¬𝑝 5 (if you are 18-25, then you are not over 25)

Acknowledgements. This work was supported by grants ANR JCJC 2019 "PRELAP" (ANR-19-CE48-0006), PEPR "iPOP" (ANR-22-PECY-0002) and ANR France 2030 "CyberINSA" (ANR-23-CMAS-0019).