A K Diongue

D Guégan

R Wolff

Exact Maximum Likelihood estimator for the BL-GARCH model under elliptical distributed innovations

Keywords: BL-GARCH process, Leverage effects, Maximum likelihood, Monte Carlo method, Volatility clustering. PACS: 65C05, 82B80, 37M10, 91B84

We are interested in the parametric class of Bilinear GARCH (BL-GARCH) models which are capable of simultaneously capturing the well known properties of financial retrun series, volatility clustering and leverage effects. Specifically, as it is often observed that the distribution of many financial time series data has heavy tails, heavier than the Normal distribution, we examine, in this paper, the BL-GARCH model in a general setting under some non-normal distributions. We also propose and implement a maximum likelihood estimation (MLE) methodology for parameter estimation. To evaluate the small-sample performance of this method for various models, a Monte Carlo study is conducted. Finally, the capability of within-sample estimation, using the S&P 500 daily returns, is also studied.

Introduction

The stylized facts of time varying variances have been long recognized in the literature. A popular and prominent tool used to describe this phenomenon is the autoregressive conditional heteroskedasticity (ARCH) model. The ARCH model, developed by [START_REF] Engle | Autoregressive conditional heteroskedasticity with estimates of the variance of the U.K. inflation[END_REF] and extended to the GARCH model by [START_REF] Bollerslev | Generalized Autoregressive Conditional Heteroskedasticity[END_REF], formulates the conditional variance of a random variable as a linear function of its past squared realizations. After the seminal work of [START_REF] Engle | Autoregressive conditional heteroskedasticity with estimates of the variance of the U.K. inflation[END_REF], a number of applied works were published to illustrate the usefulness of the GARCH models in economic and financial areas, as in [START_REF] Gouriéroux | ARCH model and financial application[END_REF] and recently [START_REF] Giraitis | Recent advances in ARCH modelling[END_REF]. Moreover, it is well established in empirical finance that the volatility of many financial assets is asymmetric. Particularly, it has been observed that stock price changes are negatively correlated with changes in volatility which implies that volatility is higher after negative shocks than after positive shocks of the same magnitude: see [START_REF] Black | Studies in stock price volatility changes[END_REF]. In the light of this empirical finding, various models with asymmetry in volatility have been proposed. Among others, they include the exponential GARCH (EGARCH) of [START_REF] Nelson | Conditional heteroskedasticity in asset returns: A new approach[END_REF], the APARCH model, proposed by [START_REF] Ding | A long memory property of stock market returns and a new model[END_REF], and the BL-GARCH model recently introduced by Storti and Vitale (2003a). The BL-GARCH in Storti and Vitale (2003a) parallel the Gaussian bilinear (BL) model of means dealt with Granger and Andersen (1978) which is used in such fields as signal, environmental studies, demography, economics and finance: see Tong (1990) for a classical review.

The present work focuses on the bilinear model in volatility (BL-GARCH), which has been investigated recently by Storti and Vitale (2003a) for its statistical properties, in that they have derived conditions for the positivity of the conditional variance and for the second order stationarity of the general model. In particluar, for the BL-GARCH(1, 1), they have provided analytical expressions for the autocorrelation function of the squared process, the unconditional fourth moment and the kurtosis coefficient, together with conditions for the existence of the higher order moments of the process. Furthermore, in order to compare the ability of the model to estimate the volatility of financial time series data with other asymmetric GARCH models, they have applied the model to the continuously compounded returns S&P 500 index. However, parameter estimation of the BL-GARCH model remains problematic. Indeed, to deal with the parameter estimation for the BL-GARCH model, Storti and Vitale (2003b) used an indirect maximum likelihood procedure based on the EM algorithm. Despite the fact that the EM algorithm has become a very popular computational method in statistics, this approach presents some limitations: slow numerical convergence, the convergence to a maximum likelihood estimator depending on a judicious choice of the starting value, and the underlying assumption of normally distributed data. These would seem to limit its applicability and also there is the issue of the non-existence of a measure of the standard errors of the estimates. In addition, while, BL-GARCH models adequately capture volatil-ity clustering and asymmetry, coupled with the auxiliary assumption of normally distributed errors, the model does not fully capture the very fat-tailed property of high frequency financial time series, often observed in exchange rates, stock returns and commodity returns, as in Baillie and Bollerslev (1991), [START_REF] Mcaleer | Automated inference and learning in modeling financial volatility[END_REF], Zivot and Wang (2005) and references therein. Actually, several alternative distributions have been proposed in the literature to model such excess kurtosis in the conditional distribution of returns better including the Student-t of [START_REF] Bollerslev | A conditional heteroscedastic time series model for speculation prices and rates of return[END_REF] and the generalized error distribution (GED) of [START_REF] Nelson | Conditional heteroskedasticity in asset returns: A new approach[END_REF]. [START_REF] Hansen | Autoregressive conditional density estimation[END_REF] introduced the use of an asymmetric Student-t type distribution to capture the skewness property well. Since then, other articles have studied different skew Student-t type distributions for financial and other applications: [START_REF] Azzalini | Distributions generated by perturbation of symmetry with emphasis on a multivariate skew t student distribution[END_REF], [START_REF] Jones | A skew extension of the t distribution with applications[END_REF] and [START_REF] Patton | On the out-of-sample importance of skewness and asymmetric dependence for asset allocation[END_REF]. All of these distributions try to account for substantial departure from the Normal distribution, specifically large kurtosis evident in the empirical distribution of the returns, more peakedness and with fatter tails than the Normal distribution. However, while large kurtosis of the returns is a well-established fact, the situation is much more obscure with regard to the symmetry of the distribution. Indeed, testing the symmetry of unconditional distributions of eight international stock market returns, Peiró (1999) concluded that, under the alternative of non-normal distributions, the symmetry of the returns cannot be rejected for most markets. Furthermore, investigating distribution-free methods, he finds that in most markets, daily financial returns are symmetric or, at least, do not present strong evidence of skewness.

The main purpose of the present study is first, following the empirical evidence of leptokurtosis in financial returns, to present an extension of the BL-GARCH with Normal distribution to other alternative non-Gaussian distributions. We therefore consider the elliptical distribution family, as it is rich enough to include the Normal, Student-t, GED and many others, [START_REF] Bingham | Semi-parametric modelling in finance: theoretical foundations[END_REF]. Secondly, given that the maximum likelihood method is widely used when dealing with parameter estimation for GARCH models, Li et al. (2002) or McNeil et al. (2005), we also investigate the problem of estimating parameters of the BL-GARCH model using MLE and examine its finite sample properties in a series of Monte Carlo simulation experiments. This is completely new, in that Storti and Vitale (2003b) only give an empirical approach for the EM method.

The paper is organized as follows. In Section 2, the BL-GARCH model is presented with some important properties concerning conditions for the conditional variance to be finite, as well as for strictly stationarity and ergodicity solutions. In Section 3, we present the MLE method under elliptical distributions, such as the Normal, Student-t and GED. We also provide the score functions as well as the Hessian matrices for these models. Section 4 studies the performance of the estimation procedure through Monte Carlo simulations. Section 5 presents the data and contains the main empirical findings along with the goodness-of-fit tests, while Section 6 provides concluding remarks.

BL-GARCH model and its specification

Over the past two decades, enormous efforts have been devoted to modelling and forecasting the volatility of stock returns and other financial time series. Seminal work in this area of research can be attributed to [START_REF] Engle | Autoregressive conditional heteroskedasticity with estimates of the variance of the U.K. inflation[END_REF], who introduced the standard autoregressive conditional heteroscedasticity (ARCH) model. Since then, further extensions have been investigated including the GARCH model proposed by [START_REF] Bollerslev | Generalized Autoregressive Conditional Heteroskedasticity[END_REF], and asymmetric GARCH models such as the EGARCH model of [START_REF] Nelson | Conditional heteroskedasticity in asset returns: A new approach[END_REF], the APARCH model of [START_REF] Ding | A long memory property of stock market returns and a new model[END_REF], and recently the BL-GARCH model, with which we deal in this paper, studied by Storti and Vitale (2003a).

The BL-GARCH model

Let S t denote the price of an asset at time t, y t = log (S t /S t-1 ) the continuously compounded return series and μ t = E (y t | Ψ t-1 ) the conditional mean given an increasing sequence of σ -fields Ψ t-1 generated by (y t-1 , y t-2 , • • •). Assume that the series of interest, y t , is given by, for t = 1, • • • , n,

y t = μ t + u t , ( 1a
)
u t = h t ε t , ( 1b
)
h 2 t = a 0 + p ∑ i=1 a i u 2 t-i + q ∑ j=1 b j h 2 t-j + r ∑ k=1 c k h t-k u t-k , ( 1c
)
where p, q, r are non-negative integers with r = min (p, q), h2 t the conditional variance of the process (u t ) t given the σ -fields Ψ t-1 and ε t ∼ i.i.d. D (0, 1) with D (.) a probability density function with mean 0 and unit variance (it will be defined in Section 3). The model (1b)-(1c) is more general than the standard GARCH model of [START_REF] Bollerslev | Generalized Autoregressive Conditional Heteroskedasticity[END_REF] in the sense that it allows innovations of different signs to have a different impact on volatility and allows larger shocks to have a larger impact on volatility than the standard GARCH model 2 .

Regularity conditions 2.2.1 Positivity conditions

Sufficient conditions for the positivity of the conditional variance h 2 t have been provided in the paper by Storti and Vitale (2003a). To this end, let

R = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ a 0 0 • • • • • • 0 0 R 1 0 • • • 0 0 . . . . . . 0 . . . . . . . . . . . . 0 • • • • • • 0 R r ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ and S = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ s 1 • • • 0 • • • 0 0 . . . 0 • • • 0 . . . . . . . . . . . . . . . 0 . . . s m-1 0 0 • • • • • • 0 s m ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 2
)
where,

for i = 1, • • • , r, R i = ⎡ ⎣ a i 1 2 c i 1 2 c i b i ⎤ ⎦ and, for i = 1, • • • , m, with m = max (p, q)- r, we set s i = a r+i , if p > q, or s i = b r+i , if p < q.
Hence, if a 0 > 0, a sufficient condition for the positivity of h 2 t is given by s i > 0 (i = 1, • • • , m) and all matrices R i (i = 1, • • • , r) being positive semi-definite. Furthermore, for the BL-GARCH(1, 1) model, the positive definiteness of the matrix R 1 is a necessary and sufficient condition for the positivity of h 2 t .

Covariance stationarity conditions

Storti and Vitale (2003a) show also that the stationarity condition of the BL-GARCH model is analogous to that stated for GARCH models by [START_REF] Bollerslev | Generalized Autoregressive Conditional Heteroskedasticity[END_REF]. Hence, the shocks h t ε t are covariance stationary and non-degenerate if

(i) var(ε t ) < ∞; (ii) a 0 > 0; and (iii) ∑ p i=1 a i + ∑ q j=1 b j < 1.
Under the assumption that the innovations are elliptically distributed with zero mean and unit variance, the first condition is obviously satisfied by construction.

The second and third conditions can either be enforced during estimation or verified afterwards.

Maximum likelihood approach

Estimation of conditional volatility models are typically performed by an MLE procedure, as in Bollerslev and Wooldridge (1992). The conditional likelihood function is given by

L (ω) = L (y 1 , • • • , y n | ω) = n ∏ t=1 g (y t , μ t (α) , h t (ω)) , (3)
where g (y t , μ t (α) , h t (ω)) denotes the conditional density function for the random variables y t with mean μ t and standard deviation h t , and ω = (α, θ ) is the parameter vector to be estimated, and where α corresponds to the set of parameters in the condtional mean assumed, in what follows, to be an ARMA(k, l) model and Storti and Vitale (2003a), we make the following theoretical assumption about the parameter space, Ω, and the true parameter vector, ω 0 = (α 0 , θ 0 ) : Assumption 1: the parameter ω 0 ∈ Ω ⊆ R k+l+p+q+r+1 is in the interior of Ω, a compact parameter space. Specifically for any vector ω ∈ Ω, assume that (1) the AR and MA polynomials have no common roots and that all their roots lie outside the unit circle;

θ = a 0 , a 1 , • • • , a p , b 1 , • • • , b q , c 1 , • • • , c r . Following
(2)

a 0 > 0, a 1 , • • • , a p ≥ 0, and b 1 , • • • , b q ≥ 0; (3) c 2 i < 4a i b i , for i = 1, • • • , r; and (4)
p ∑ i=1 a i + q ∑ i=1 b i < 1.
Thus, estimation proceeds by maximising, under the Assumption 1, L (ω) = log (L (ω)), where pre-sample values of h 2 t are set to the unconditional sample variance.

Since it may be expected that excess kurtosis and skewness displayed by the residuals of conditional heteroscedasticity models will be reduced when a more appropriate distribution is used, we consider in this study the three most typical elliptical normalized distributions that have been applied so far: the Normal, Student-t and GED distributions.

Normal distribution

The Normal distribution is the most widely used when estimating GARCH models.

If we assume that the innovations (ε t ) t∈Z have a conditional Gaussian distribution then the conditional log-likelihood function associated to y t , u t = y t -μ t , is given by

L (ω) = - n 2 log (2π) - 1 2 n ∑ t=1 log h 2 t + u 2 t h 2 t , ( 4
)
where n is the number of observations. To obtain an analytical or numerical solution for the MLE, we need the first-order derivative and to solve the equation ∂ L (ω) /∂ ω = 0. Taking the differential of L (ω) in ( 4) with respect to the full set of parameter ω yields

∂ L (ω) ∂ ω = n ∑ t=1 u t h 2 t ∂ μ t ∂ ω + 1 2 n ∑ t=1 1 h 2 t u 2 t h 2 t -1 ∂ h 2 t ∂ ω . ( 5
)
The Hessian matrix is given by

∂ 2 L (ω) ∂ ω∂ ω = - n ∑ t=1 1 h 2 t ∂ μ t ∂ ω ∂ μ t ∂ ω - n ∑ t=1 u t h 4 t ∂ μ t ∂ ω ∂ h 2 t ∂ ω - n ∑ t=1 u t h 4 t ∂ μ t ∂ ω ∂ h 2 t ∂ ω + n ∑ t=1 1 h 4 t 1 2 - u 2 t h 2 t ∂ h 2 t ∂ ω ∂ h 2 t ∂ ω . ( 6
)

Student-t distribution

Now, if we assume that the innovations (ε t ) t∈Z have a conditional Student-t distribution with ν degrees of freedom, then the MLE estimator ωn maximises the log-likelihood function L (ω) given by

L(ω) = n log Γ ν + 1 2 -log Γ ν 2 - 1 2 log π (ν -2) - 1 2 n ∑ t=1 log(h 2 t ) + (ν + 1) log 1 + u 2 t h 2 t (ν -2) , ( 7
)
where 2 < ν ≤ ∞ and Γ is the Euler gamma function defined by Γ (x) = ∞ 0 t x-1 e -t dt. When ν → ∞, we have the Normal distribution, so that the smaller the value of ν the fatter the tails. The score function is given by

∂ L (ω) ∂ ω = n ∑ t=1 ν + 1 ν -2 u t h 2 t 1 + u 2 t h 2 t (ν -2) -1 ∂ μ t ∂ ω + 1 2 n ∑ t=1 ν + 1 ν -2 u 2 t h 2 t 1 + u 2 t h 2 t (ν -2) -1 -1 1 h 2 t ∂ h 2 t ∂ ω , ( 8
)
and the Hessian matrix is given by

∂ 2 L (ω) ∂ ω∂ ω = ν + 1 ν -2 n ∑ t=1 1 + u 2 t (ν -2) h 2 t -1 u t h 2 t 2 ν -2 1 + u 2 t (ν -2) h 2 t -1 u t h 2 t -1 ∂ μ t ∂ ω ∂ μ t ∂ ω + ν + 1 ν -2 n ∑ t=1 1 + u 2 t (ν -2) h 2 t -1 u t h 4 t 1 ν -2 1 + u 2 t (ν -2) h 2 t -1 u 2 t h 2 t -1 ∂ μ t ∂ ω ∂ h 2 t ∂ ω + ν + 1 ν -2 n ∑ t=1 1 + u 2 t (ν -2) h 2 t -1 u t h 4 t 1 ν -2 1 + u 2 t (ν -2) h 2 t -1 u 2 t h 2 t -1 ∂ h 2 t ∂ ω ∂ μ t ∂ ω + 1 2 n ∑ t=1 1 h 4 t 1 + (ν + 1) u 2 t (ν -2) h 2 t 1 + u 2 t (ν -2) h 2 t -1 . u 2 t (ν -2) h 2 t 1 + u 2 t (ν -2) h 2 t -1 -2 ∂ h 2 t ∂ ω ∂ h 2 t ∂ ω . ( 9
)

GED distribution

Knowing that skewness and kurtosis are important in financial applications, [START_REF] Nelson | Conditional heteroskedasticity in asset returns: A new approach[END_REF] suggested to consider the family of GEDs. The probability density function of a normalized GED random variable is given by

f (x) = ν2 -(1+ 1 ν ) λ ν Γ 1 ν e -1 2 x λ ν ν , -∞ < x < ∞, (10)
with λ ν = 2 -2/ν Γ (1/ν) /Γ (3/ν) and 0 < ν < ∞ is the tail-thickness parameter.

The GED includes the Gaussian distribution (ν = 2) as a special case, along with many other distributions, some more fat-tailed than the Gaussian (e.g., the double exponential distribution corresponding to ν = 1) and some more thin-tailed (e.g., the Uniform distribution on the interval -√ 3, √ 3 when ν → ∞). The GED loglikelihood function is given by

L(ω) = n log ν λ ν -1 + 1 ν log (2) -log Γ 1 ν - 1 2 n ∑ t=1 log h 2 t + h -ν t u t λ ν ν . ( 11
)
In this case, the score function and the Hessian matrix, respectively, are given by

∂ L (ω) ∂ ω = 1 2 n ∑ t=1 ν |λ ν | u t h t ν 1 u t ∂ μ t ∂ ω + 1 2 n ∑ t=1 1 h 2 t 1 2 l |λ ν | ν u t h t ν -1 ∂ h 2 t ∂ ω , (12)
and

∂ 2 L (ω) ∂ ω∂ ω = - ν |λ ν | ν ν -3 2 n ∑ t=1 1 h 2 t u 2 t h 2 t ν 2 -1 ∂ μ t ∂ ω ∂ μ t ∂ ω - 1 2 ν |λ ν | ν n ∑ t=1 u t h 4 t u 2 t h 2 t ν 2 -1 1 + ν 2 -1 u 2 t h 2 t - 1 2 ν |λ ν | ν n ∑ t=1 u t h 4 t u 2 t h 2 t ν 2 -1 1 + ν 2 -1 u 2 t h 2 t - 1 2 n ∑ t=1 1 4 ν (ν + 2) |λ ν | ν u t h t ν -1 1 h 4 t ∂ h 2 t ∂ ω ∂ h 2 t ∂ ω . ( 13
)
Equations ( 5) through [START_REF] Mcaleer | Automated inference and learning in modeling financial volatility[END_REF] require the computation of ∂ h 2 t /∂ ω. This is given by

∂ h 2 t ∂ ω = 1, u 2 t-1 , • • • , u 2 t-p , h 2 t-1 , • • • , h 2 t-q , h t-1 u t-1 , • • • , h t-r u t-r +2 p ∑ i=1 a i u t-1 ∂ u t-i ∂ ω + q ∑ j=1 b j ∂ h 2 t-j ∂ ω + r ∑ k=1 c k h t-k ∂ u t-k ∂ ω + 1 2 u t-k h t-k ∂ h 2 t-k ∂ ω . ( 14
)
The use and analysis of the MLE method for the estimation problem is classical. A main attraction is the general feature that maximum likelihood estimators achieve optimal accuracy, in that they are asymptotically consistent (in data length n), and achieve the Cramér-Rao lower bound on estimate variability. Despite these advantages, an important obstacle to employing the method is the difficulty of computing a value ωMLE that satisfies Assumption 1. In the next section, the practical applicability and small sample performance of the MLE procedure for BL-GARCH processes are studied by Monte Carlo simulations.

Monte Carlo experiments

To our knowledge, no results exist on the properties of these estimators when we observe a finite segment of (1a)-(1c). Thus, we have designed and executed a Monte Carlo experiment using the different distributions described in the previous section as data generating processes, with the aim of analyzing the sampling properties of the exact MLE estimators of the parameter vector ω for the BL-GARCH model. Through the Monte Carlo experiment, the model considered for u t = y t -μ t is a BL-GARCH(1, 1) given by

u t = h t ε t , for t = 1, • • • , n, ( 15a
)
h 2 t = a 0 + a 1 u 2 t-1 + b 1 h 2 t-1 + c 1 h t-1 u t-1 , (15b)
with ε t is a standard Normal, Student-t or GED random variable and n = 500, 1000, 2000 and 3000. Two cases are studied in the simulation experiments. In the first case, the conditional mean, μ t , is taken equal to zero while in the second one we assume that it follows an AR(1) model 3 . The data generating processes are summarized in Table 1 with the first three lines corresponding to the case μ t = 0 and the last line μ t = α 0 + α 1 y t-1 . The data generating processes, correponding to the case μ t = 0, are as in Storti and Vitale (2003a). Throughout the simulations, we consider a Student-t with five degrees of freedom, the same as used in [START_REF] Tsay | Analysis of financial time series[END_REF]. Thus, the first four moments of the conditional density exist. For the GED distribution, we assume that the tail-thickness parameter is equal to three. 3 and4 present the estimated results from a non-Gaussian BL-GARCH model. We read, from these tables, that the averages of the parameter estimates are close to the true values under each of the underlying non-Normal error terms. The MAE as well as the RMSE are quite small and decrease when the sample sizes increase. Finally, Table 5 summarizes the results from the AR(1)-BL-GARCH(1, 1). Results reveal that parameter estimates are satisfactory in that the MAE and also RMSE are small. We can also remark that, in general, the estimators of the autoregressive models seem not to be affected by, the presence of the BL-GARCH errors. In addition, the method seems applicable, even if the sample size is less than 100, due to the fact that, in general, the true values are contained with ±2 standard deviations of MLE's estimates: see Figures 1234.

Normal

Empirical study

The daily continuously compounded returns of the S&P 500 stock market index are used for the empirical study in this paper to gauge the effectiveness of the BL-GARCH-type model with Normal, Student-t and GED innovations. In particular, we analyze the period from March 01, 1999 through January 31, 2001, which yields n = 487 daily observations, excluding public holidays. The sample closely corresponds to the data used by Storti and Vitale (2003b). Table 6 gives the summary statistics of the S&P 500 log returns for the full sample. The mean and the standard deviation are quite small, while the estimated measure of skewness is significantly positive, indicating that the S&P 500 has non-symmetric returns. The kurtosis is a little higher than that of a Normal distribution which is 3, suggesting that fat-tailed distributions could better describe the unconditional distribution of the data. The results of the non-Normality test agree with prior literature using financial data, that is, a leptokurtic distribution is found for the S&P 500 log return data. The Box-Pierce Q-tests of up to twenty-fourth order serial correlation for the levels and squares of the mean-corrected S&P 500 log returns were performed. Q (24) and Q 2 (24) are significant for both the return and squared return series. The diagnostics suggest that a GARCH-class model would be appropriate, along with an error distribution that allows for greater kurtosis than the Normal distribution. The Jarque Bera test critical value at significance level of 5% is 5.85423. Q and Q 2 are the Box Pierce statistics for the levels and squared of the S&P 500 log returns respectively, using 24 lags. The critical value at level of 5% is 36.4150.

Figure 5 gives, the time plot of the data while Figure 6 shows that the returns distribution also exhibits fat-tails confirming the results in Table 6. 7, are the log-likelihood function at its maximum, the Akaike Information Criterion (AIC) and the Bayesian Information Crietrion (BIC). We further report the values of the Box-Pierce (Q) statistics for the standardized and squared standardized residuals with in parentheses, the corresponding p-values as a check of the empirical validity of the models. It is clear from Table 7 that estimate â1 and b1 in the GARCH(1, 1) are significant at the 5% level with the volatility coefficient greater in magnitude. Hence the hypothesis of constant variance can be rejected, at least within sample. Furthermore, the stationarity condition is satisfied for the three distributions, as â1 + b1 < 1 at the maximum of the respective log-likelihood functions.

From results obtained on estimation of the BL-GARCH(1, 1) model, we remark, once again, that the innovation and volatility spillovers are significantly different from zero. Further, the estimated asymmetric volatility response ( ĉ1 ) is negative and significant for all models confirming the usual expectation in stock markets where downward movements (falling returns) are followed by higher volatility than upward movements (increasing returns). In all cases, the tails parameter estimate is strongly significant far two for the GED distribution while it is very large through insignificant for the Student-t distribution but 1/ν is also different to zero. However, the results confirm the empirical findings by Storti and Vitale (2003a), in that the kurtosis strongly depends to the leverage-effect response parameter.

The results for Q-statistics shown in Table 7 are not significant up to order 12 and also order 24, which indicates that the BL-GARCH(1, 1) as well as the GARCH(1, 1) process are appropriate to model the conditional variance of the S&P 500 logreturns. However, the goodness-of-fit statistics as well as the residuals diagnostics indicate that the BL-GARCH performs better in describing the conditional variance of the S&P 500 returns. Moreover, the possible usefulness of using fat-tailed innovations for the BL-GARCH model seems to be confirmed by the log-likelihood values and the AIC.

Conclusion

This study obtains exact maximum likelihood estimates of a BL-GARCH process with conditionally elliptical distributions. The small-sample properties indicate that the approach can yield asymptotically efficient estimates. In addition, these results strongly suggest that the maximum likelihood estimation inference procedure can be used to estimate the parameters of the BL-GARCH model, even in samples as small as 100 observations. Further, as the simulation experiments show, one advantage of the maximum likelihood estimator procedure, proposed in this paper, compared to the method used by Storti and Vitale (2003b), is that it could simultaneously estimate the parameters of the BL-GARCH when the conditional mean is assumed non-constant. Further, the empirical results reveal that the BL-GARCHt(1, 1), i.e a BL-GARCH model with conditional errors that are t-distributed, fits the data best. This table provides the estimated coefficients, standard errors for the conditional standard deviation equation for the S&P 500 log returns index market. â0 is the constant in the conditional standard deviation equation, â1 is the ARCH coefficient, b1 is the GARCH coefficient, ĉ1 is the leverage effect, ν is the degrees of freedom. Log-lik is the maximized log likelihood. AIC is the Akaike Information Citerion and BIC the Bayesian Information Criterion. Q and Q 2 are the Box Pierce statistics for the standardized and squared standardized residuals respectively, using 12 and 24 lags with p-values in square brackets. The critical values at significant level of 5% are 21.026069 and 36.415028 respectively.

Fig. 1 . 1 Fig. 2 . 1 Fig. 3 . 1 Fig. 4 .

 1121314 Fig. 1. Boxplot of the parameter a 0 of Model 1 under Normal, Student-t and GED distributions (n = 1000)

Table 1

 1 This operation is performed in order to avoid dependence on initial values. The calculations were carried out in Matlab on a Pentium IV CPU 3.00 GHz computer. Inspection of Table2, corresponding to the Normal case, reveals that, for both sample sizes, the averages obtained from the exact MLE are close to the true parameter values. The corresponding MAE and RMSE are very small indicating that estimators are asymptotically unbiased and consistent, respectively. Tables

	Data generating processes (DGPs)
	DGP	α 0	α 1	a 0	a 1	b 1	c 1
	Model 1	0	0 0.01 0.09 0.9 0.15
	Model 2	0	0 0.05 0.05 0.9 0.25
	Model 3	0	0	0.2 0.05 0.75 0.35
	Model 4 0.01 0.2 0.01 0.09 0.9 0.15

Tables 2-4 list the Monte Carlo mean, mean absolute error (MAE) and root mean square error (RMSE) for the parameter vector ω across M = 1000 Monte Carlo simulations. The simulation algorithm generates n + 500 observations for each series, saving only the last n.

Table 2 .

 2 Estimator parameters for the centered Gaussian BL-GARCH model defined by (15a)-(15b)

	â0 â1 b1 ĉ1	Model n Mean MAE RMSE Mean MAE RMSE Mean MAE RMSE Mean MAE RMSE	1 100 0.09565 0.08864 0.26478 0.09236 0.07241 0.09654 0.76343 0.18036 0.29031 0.15901 0.09549 0.13342	300 0.01770 0.01032 0.02273 0.08398 0.03705 0.04610 0.89119 0.04473 0.07083 0.15306 0.03807 0.04939	1000 0.01150 0.00308 0.00422 0.08943 0.01529 0.01953 0.89718 0.01557 0.02021 0.15252 0.01951 0.03658	3000 0.01047 0.00156 0.00205 0.09022 0.00875 0.01090 0.89883 0.00874 0.01089 0.15129 0.01124 0.01463	2 100 0.09694 0.06658 0.12746 0.06368 0.05110 0.08043 0.82400 0.12195 0.21326 0.27502 0.09176 0.12345	300 0.05385 0.01643 0.02661 0.04573 0.02613 0.03235 0.89923 0.03529 0.05210 0.25481 0.04424 0.05712	1000 0.05141 0.00741 0.00972 0.04756 0.01312 0.01637 0.90081 0.01644 0.02073 0.25072 0.02318 0.02904	3000 0.05010 0.00369 0.00469 0.04905 0.00719 0.00884 0.90071 0.00870 0.01073 0.24942 0.01220 0.01540	3 100 0.22124 0.09400 0.12953 0.07752 0.04822 0.07592 0.70142 0.12867 0.18151 0.36293 0.10454 0.13423	300 0.20409 0.04860 0.06385 0.05699 0.02454 0.03419 0.73837 0.06370 0.08437 0.35189 0.06176 0.07852	1000 0.20094 0.02428 0.03056 0.05201 0.01420 0.01823 0.74715 0.03134 0.03988 0.35259 0.03335 0.04184	3000 0.20013 0.01390 0.01751 0.05023 0.00884 0.01096 0.74947 0.01779 0.02252 0.35015 0.01918 0.02407	This table summarizes the estimates coefficients from the BL-GARCH(1,1) model with the true value set of parameters {a 0 ,a 1 ,b 1 ,c 1 } = {0.01,0.09,0.9,0.15}, {0.05,0.05,0.9,0.25} and	{0.2,0.05,0.75,0.35} . MAE -Mean Absolute Errors, RMSE -root mean square error for the Gaussian BL-GARCH(1,1). Monte Carlo simulations are computed with 1000 replications. Each repli-	cation gives a sample size n = 100, 200, 1000 and 3000 of observations.

Table 3 .

 3 Estimator parameters for the centered Student-t BL-GARCH model defined by (15a)-(15b)

	â0 â1 b1 ĉ1	Model n Mean MAE RMSE Mean MAE RMSE Mean MAE RMSE Mean MAE RMSE	1 100 0.10610 0.09896 0.27096 0.10400 0.08644 0.012400 0.73671 0.20345 0.31462 0.14972 0.13008 0.18310	300 0.02153 0.01403 0.04723 0.08986 0.04244 0.05539 0.88264 0.05100 0.09447 0.15937 0.04995 0.06570	1000 0.0119 0.00347 0.00473 0.09018 0.0189 0.0242 0.8961 0.0185 0.0239 0.1522 0.0256 0.0321	3000 0.0104 0.00162 0.00214 0.08951 0.0103 0.0129 0.8994 0.0095 0.0119 0.1493 0.0145 0.0184	2 100 0.14528 0.11232 0.18852 0.07323 0.05744 0.09183 0.76433 0.17098 0.27763 0.28221 0.13470 0.18218	300 0.07061 0.03569 0.08397 0.05241 0.03002 0.04431 0.87763 0.06095 0.12948 0.26215 0.06112 0.07734	1000 0.05164 0.00949 0.01263 0.04824 0.01267 0.01606 0.90045 0.01684 0.02198 0.25186 0.03032 0.03909	3000 0.05058 0.00525 0.00624 0.04958 0.00727 0.00906 0.90005 0.00959 0.01209 0.25087 0.01801 0.02241	3 100 0.24646 0.13090 0.25177 0.08733 0.06070 0.10274 0.67189 0.16486 0.23289 0.35433 0.13832 0.17808	300 0.21268 0.07507 0.09401 0.06466 0.03185 0.04799 0.71045 0.10450 0.13967 0.35592 0.07987 0.10375	1000 0.2018 0.03793 0.04975 0.05065 0.01584 0.02003 0.74725 0.04575 0.05793 0.34400 0.04590 0.05716	3000 0.20206 0.02003 0.02479 0.04994 0.00940 0.01191 0.74910 0.02282 0.02897 0.35035 0.02673 0.03202	This table summarizes the estimates coefficients from the BL-GARCH(1,1) model with the true value set of parameters {a 0 ,a 1 ,b 1 ,c 1 } = {0.01,0.09,0.9,0.15}, {0.05,0.05,0.9,0.25} and	{0.2,0.05,0.75,0.35} . MAE -Mean Absolute Errors, RMSE -root mean square error for the Student-t BL-GARCH(1,1). Monte Carlo simulations are computed with 1000 replications. Each repli-	cation gives a sample size n = 100, 200, 1000 and 3000 of observations.

Table 4 .

 4 Estimator parameters for the centered GED BL-GARCH model defined by (15a)-(15b)

	â0 â1 b1 ĉ1	Model n Mean MAE RMSE Mean MAE RMSE Mean MAE RMSE Mean MAE RMSE	1 100 0.07398 0.06648 0.14326 0.10136 0.07136 0.08697 0.77613 0.15628 0.24907 0.16339 0.08021 0.10712	300 0.01609 0.00806 0.01114 0.08845 0.02841 0.03609 0.89029 0.03428 0.04430 0.15484 0.03144 0.03870	1000 0.01105 0.00261 0.00351 0.08917 0.01448 0.01832 0.89884 0.01432 0.01823 0.15094 0.01513 0.01919	3000 0.01025 0.00132 0.00181 0.09067 0.00778 0.01027 0.89866 0.00778 0.01026 0.15070 0.00894 0.01183	2 100 0.08216 0.04731 0.09275 0.06423 0.05041 0.07831 0.84761 0.09293 0.16582 0.27555 0.07683 0.10431	300 0.05405 0.01372 0.01935 0.04715 0.02526 0.03173 0.89845 0.03165 0.04357 0.25623 0.03677 0.04708	1000 0.05077 0.00567 0.00738 0.04875 0.01225 0.01547 0.90018 0.01433 0.01819 0.25067 0.01906 0.02358	3000 0.05037 0.00312 0.00393 0.04956 0.00684 0.00866 0.90008 0.00782 0.00976 0.25066 0.01022 0.01287	3 100 0.21931 0.08345 0.11359 0.07424 0.04509 0.07298 0.70846 0.11736 0.16225 0.36299 0.09512 0.12272	300 0.21268 0.07507 0.09401 0.06466 0.03185 0.04799 0.71045 0.10450 0.13967 0.35592 0.07987 0.10375	1000 0.20262 0.02017 0.02578 0.05130 0.01263 0.01631 0.74550 0.02690 0.03466 0.35097 0.02695 0.03392	3000 0.20020 0.01124 0.01393 0.05082 0.00833 0.01040 0.74919 0.01494 0.01874 0.35114 0.01548 0.01958	This table summarizes the estimates coefficients from the BL-GARCH(1,1) model with the true value set of parameters {a 0 ,a 1 ,b 1 ,c 1 } = {0.01,0.09,0.9,0.15}, {0.05,0.05,0.9,0.25} and	{0.2,0.05,0.75,0.35} . MAE -Mean Absolute Errors, RMSE -root mean square error for the GED BL-GARCH(1,1). Monte Carlo simulations are computed with 1000 replications. Each replica-	tion gives a sample size n = 500, 1000, 2000 and 3000 of observations.

Table 5

 5 Estimator parameters for the centered AR-BL-GARCH model defined by (15a)-(15b)

	Distribution	n	α0	α1	â0	â1	b1	ĉ1
	Normal	100	0.0257	0.18447	0.10767	0.09233	0.78056	0.16659
			(0.06246)	(0.08631)	(0.08060)	(0.07867)	(0.17161)	(0.09643)
			[0.09254]	[0.08737]	[0.24972]	[0.10633]	[0.28043]	[0.13219]
		300	0.01897	0.19377	0.01874	0.0838	0.88978	0.15629
			(0.03242)	(0.04678)	(0.01141)	(0.03729)	(0.04555)	(0.03919)
			[0.04338]	[0.05949]	[0.02728]	[0.04697]	[0.07751]	[0.05055]
		1000	0.01103	0.19871	0.01142	0.08847	0.89841	0.15103
			(0.01624)	(0.02544)	(0.00322)	(0.01555)	(0.01603)	(0.01932)
			[0.02069]	[0.03167]	[0.00441]	[0.01995]	[0.02064]	[0.02422]
		3000	0.01117	0.19905	0.01042	0.08998	0.89941	0.15086
			(0.00911)	(0.01477)	(0.00156)	(0.0088)	(0.00871)	(0.01129)
			[0.01189]	[0.0187]	[0.00209]	[0.01134]	[0.01117]	[0.01444]
	Student-t	100	0.0209	0.18876	0.08852	0.09884	0.7457	0.15474
			(0.05193)	(0.08171)	(0.08201)	(0.08648)	(0.20145)	(0.13304)
			[0.07803]	[0.10338]	[0.23158]	[0.11957]	[0.32115]	[0.17796]
		300	0.01472	0.19351	0.02124	0.08836	0.87681	0.15205
			(0.02583)	(0.04575)	(0.01389)	(0.04023)	(0.05523)	(0.05247)
			[0.03494]	[0.05758]	[0.03351]	[0.05389]	[0.10593]	[0.06717]
		1000	0.01123	0.19838	0.01206	0.09095	0.89505	0.15105
			(0.01378)	(0.02474)	(0.0036)	(0.01713)	(0.01776)	(0.02524)
			[0.01737]	[0.031]	[0.00523]	[0.02215]	[0.0227]	[0.03156]
		3000	0.01053	0.19941	0.01054	0.09075	0.89861	0.15111
			(0.00794)	(0.01389)	(0.0018)	(0.00992)	(0.00957)	(0.0151)
			[0.01013]	[0.01748]	[0.00278]	[0.01286]	[0.01323]	[0.01999]
	GED	100	0.01715	0.18247	0.09296	0.09692	0.77744	0.16364
			(0.06385)	(0.08462)	(0.08650)	(0.07818)	(0.17089)	(0.08186)
			[0.09838]	[0.10735]	[0.10735]	[0.10405]	[0.28519]	[0.11175]
		300	0.01680	0.19660	0.01619	0.08756	0.89222	0.15398
			(0.03025)	(0.04396)	(0.00882)	(0.03510)	(0.03925)	(0.0323)
			[0.04067]	[0.05591]	[0.02660]	[0.04534]	[0.06616]	[0.04133]
		1000	0.01253	0.19973	0.01098	0.0914	0.89718	0.15228
			(0.0163)	(0.02493)	(0.0028)	(0.01539)	(0.01545)	(0.01624)
			[0.021]	[0.03152]	[0.004]	[0.02045]	[0.02062]	[0.02077]
		3000	0.01059	0.19945	0.01026	0.09187	0.89782	0.15056
			(0.00949)	(0.01377)	(0.00159)	(0.00928)	(0.0914)	(0.01029)
			[0.01288]	[0.01753]	[0.00232]	[0.01533]	[0.01465]	[0.01379]
	This table summarizes the estimates coefficients from the AR(1)-BL-GARCH(1,1) model with the true value set of parame-
	ters {α 0 ,α 1 ,a 0 ,a 1 ,b 1 ,c 1 } = {0.01,0.2,0.01,0.09,0.9, 0.15}. MAE -Mean Absolute Errors, RMSE -root mean square error
	for the AR(1)-BL-GARCH(1,1). Monte Carlo simulations are computed with 1000 replications. Each replication gives a
	sample size n = 500, 1000, 2000 and 3000 of observations.

Table 6

 6 Statistics of daily log returns of the S&P 500 stock market index.

	Number of observations 487	Skewness	0.03708
	Mean	0.0002055	Kurtosis	4.4643
	Standard deviation	0.01283	Jarque Bera test 43.3524
	Minimum	-0.06004	Q (24)	39.2169
	Maximum	0.04888	Q 2 (24)	45.7208

Table 7

 7 presents the results of maximum likelihood estimation of the parameters of the model under candidate the various assumptions on innovations, with the standard errors computed using the Hessian matrix4 . The parameter estimates from Fig.6. Non-parametric density of S&P 500 daily returns and probability density function of the normal distribution, in dotted line the GARCH model with Normal, Student-t and GED errors are also provided in this table. In order to compare objectively, various goodness-of-fit statistics are used. The diagnostics, summarized in Table

		0.04
		0.02
	Daily returns	-0.02 0
		-0.04
		-0.06
		0	50	100	150	200	250	300	350	400	450	500
							Time
			Fig. 5. S&P 500 daily returns 03/01/99-01/31/01

Table 7

 7 Conditional variance model estimates for the S&P 500 return

		GARCH model		BL-GARCH model
	Parameters	Normal	Student-t	GED	Normal	Student-t	GED
	â0	0.00000715	0.000006361	0.000006108	0.000011394	0.000009243	0.00001057
		(0.28104 10 -10 )	(0.20586 10 -10 ) (0.20144 10 -10 )	(0.14775 10 -10 )	(0.13496 10 -10 )	(0.15424 10 -10 )
	â1	0.0568184	0.0504452	0.0499355	0.060119	0.0513906	0.0559965
		(0.0007625)	(0.0006168)	(0.0006111)	(0.0006226)	(0.0005249)	(0.0006235)
	b1	0.900428	0.911250	0.913097	0.880531	0.900687	0.888784
		(0.0026979)	(0.0017480)	(0.0018101)	(0.0013771)	(0.0011259)	(0.0013873)
	ĉ1	-	-	-	-0.271323	-0.249673	-0.261943
		(-)	(-)	(-)	(0.0027800)	(0.0030085)	(0.0031051)
	ν	-	8.676471	1.514841	-	14.943269	1.741412
		(-)	(9.929491)	(0.017712)	(-)	(67.446372)	(0.024844)
				Goodness-of-fit statistics
	Log-lik	1435.91706	1441.94452	1441.04650	1456.47965	1458.63396	1457.65676
	AIC	-2865.83441	-2875.88904	-2874.09300	-2904.95930	-2907.26792	-2905.31353
				Diagnostics
	Q (12)	15.892571	16.018537	16.010225	15.876516	15.739415	15.833319
		(0.196206)	(0.190388)	(0.190768)	(0.196958)	(0.203465)	(0.198990)
	Q (24)	32.082112	32.206641	32.159427	32.314595	32.501633	32.417392
		(0.124969)	(0.121951)	(0.123088)	(0.119381)	(0.115033)	(0.116975)
	Q 2 (12)	4.682633	4.589348	4.566436	4.477202	4.455119	4.485078
		(0.967748)	(0.970309)	(0.970918)	(0.973208)	(0.973756)	(0.973011)
	Q 2 (24)	19.326186	19.483667	19.377018	22.177684	22.542585	22.475412
		(0.734378)	(0.725728)	(0.731595)	(0.568662)	(0.546916)	(0.550914)

Changes in stock prices tend to be negatively correlated related to changes in volatility,[START_REF] Black | Studies in stock price volatility changes[END_REF]. Specifically, lagged unexpected declines in prices (bad news) increase current volatility more than to lagged unexpected increases in prices (good news). This asymmetry in stock returns volatility has been termed the "leverage-effect".

Other simulation results are available upon request.

The time for convergence of the algorithm under specific distribution is approximately

seconds.

This research was supported by an Australian Research Council Discovery Project grant (DP0666266).


domain: math.math-ds math sdv

text:
Jairo Gomes
email: jairo.gomes@ifmt.edu.br

Izabel Da A Silva

Cristina Rodrigues

Mostafa Adimy
email: mostafa.adimy@inria.fr

Paulo Fernando

Arruda Mancera
email: paulo.mancera@unesp.br

Izabel Cristina Rodrigues

Paulo Fernando De

J G Da Silva

I C Da Silva

P F Mancera

The effect of lenvatinib and pembrolizumab on thyroid cancer refractory to iodine 131I simulated by mathematical modeling

Keywords: Thyroid tumor, Targeted therapy, Immunotherapy, Immune system, Mathematical model, Asymptotic stability, Lyapunov function

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

Introduction

Differentiated thyroid cancer (DTC) is usually treated with surgeries such as lobectomy, thyroidectomy, and, when necessary, radioactive iodine (RAI) therapy based on 131 I. [START_REF] Haugen | American Thyroid Association management guidelines for adult patients with thyroid nodules and differentiated thyroid cancer[END_REF] However, after malignant thyroid cells lose their ability to capture and concentrate iodine, other therapeutic options such as targeted therapies and immunotherapies should be applied. [START_REF] Tirrò | Molecular alterations in thyroid cancer: from bench to clinical practice[END_REF] According to Ref. 3, clinical and pathological evidence suggest that the increasingly poor differentiation of DTC until becoming anaplastic cancer is a natural process in the development of malignancy. In addition, papillary thyroid cancer (PTC) carrying mutations in the BRAF and RAS genes and RET-PTC rearrangements are more likely to progress through stages of increasingly poor differentiation. Those genetic changes can activate the mitogen-activated protein kinase (MAPK) signaling pathway, which is responsible for regulating various cellular activities. For that reason, those mutations have been explored as therapeutic targets in cancer treatment. [START_REF] Nikiforov | Thyroid carcinoma: molecular pathways and therapeutic targets[END_REF]5 An example of drug used to treat DTC due to a loss of sensitivity to iodine is lenvatinib a . According to Ref. 2, lenvatinib is a targeted therapeutic agent that acts as an inhibitor of multiple tyrosine kinases and primarily aims to reduce the proliferation of tumor cells by blocking vascular endothelial growth factors (VEGF).

In a randomized, double-blind, phase 3 study, lenvatinib was assessed by Ref. 6. The study had a total of 392 participants, 131 of whom received the placebo. The screening revealed several important results, including that responses rate to treatment of 64.8% and 1.5% using lenvatinib and placebo, respectively; the progression of cancer in 35.6% of patients receiving lenvatinib versus 83.2% in those taking placebo during the follow-up period, being approximately 17 months in both cases; and an average tumor progression-free survival of 18.3 months in the group using lenvatinib and 3.6 months in those using placebo. Some of the side effects observed in the study were hypertension, fatigue, and diarrhea. [START_REF] Schlumberger | Lenvatinib versus placebo in radioiodine-refractory thyroid cancer[END_REF] Considering the important role of the immune system in the development of cancer, Ref. 7 examined the infiltration of regulatory T cells (Tregs) and natural killer (NK) cells in thyroid gland of patients with PTC. Among the results, some effector T cells were strongly controlled by Tregs cells that inhibited their proliferation, as well as prevented attacks on the tumor by the immune system. The infiltration of Tregs was positively associated with advanced stages of tumor progression, whereas NK cells were negatively correlated with those stages. As for lenvatinib, Ref. 8 observed that in cases of murine melanoma and kidney cancer, the infiltration of NK cells into the tumor increased as a result of targeted therapy. Thus, the drug also caused immune modulation in tumors, thereby making it a promising complementary agent to be used with immunotherapies in order to increase tumor infiltration and NK cell activation.

The immune regulatory molecule known as programmed cell death ligand 1

a Medicine developed and marketed by Eisai Co. Ltd. under the label Lenvima.

(PD-L1), correlated with PTC involving a greater risk of metastasis or shorter disease-free survival, is indicated as a prognostic marker of cancer under treatment with immunotherapy. 9, 10 PD-L1 blocks the immune response of activated T cells by binding to the PD-1 receptor, an immune checkpoint expressed in T cells, to modulate their activation or inhibition. [START_REF] Tirrò | Molecular alterations in thyroid cancer: from bench to clinical practice[END_REF][START_REF] Du Rusquec | Clinical utility of pembrolizumab in the management of advanced solid tumors: an evidencebased review on the emerging new data[END_REF] Pembrolizumab b is a monoclonal antibody of human immunoglobulin G4 able to prevent the binding of PD-1 to PD-L1 and thus restore the antitumor immune response of T cells. The antitumor response observed by using targeted therapies in PTC that is refractory to RAI (RAIR) is usually limited, in the sense that the tumor always regrows at some point during or after treatment. [START_REF] Mehnert | Safety and antitumor activity of the anti-PD-1 antibody pembrolizumab in patients with advanced, PD-L1-positive papillary or follicular thyroid cancer[END_REF] Several studies have addressed the combination of targeted therapies and immunotherapy as a treatment for RAIR-DTC. For example, ongoing clinical trials currently evaluate the combination of lenvatinib and pembrolizumab in patients with DTC (i.e., NCT02973997 and NCT02628067) 2, 12-15 and other cancers. [START_REF] Kato | Lenvatinib plus anti-PD-1 antibody combination treatment activates CD8 + T cells through reduction of tumor-associated macrophage and activation of the interferon pathway[END_REF][START_REF] Makker | Lenvatinib plus pembrolizumab in patients with advanced endometrial cancer: an interim analysis of a multicentre, open-label, single-arm, phase 2 trial[END_REF] Mathematical modeling has been widely used to study tumor dynamics and cancer treatments. For instance, considering three coupled differential equations and variables at time t denoting the tumor cell population, the total level of NK cells, and the effectiveness of tumor-specific CD8 + T cell, Ref. 18 modeled tumorimmune interactions to clarify the process of immune-mediated tumor rejection. In another study, Ref. 19, in addition to including tumor cells, host cells, and immune cells, the authors simulated drug interactions in order to identify the best protocols for chemotherapy. Immunotherapy has also been studied with mathematical models. Ref. 20 applied a system of ordinary differential equations (ODE) and the theory of optimal control to identify the best time of administration and dosage for dendritic cell vaccines designed to reduce or eliminate the tumor mass. In Ref. 21, the authors addressed the possibility of combining immunotherapy and radiotherapy based on a tumor growth model and examined the use of radiotherapy with inhibitors of the PD1-PDL1 axis and/or the CTLA4 pathway to describe how the antitumor immune response occurs. Ref. 22 present a review about mathematical models of thyroid cancer, one that discusses different models used and aspects investigated to date. In Ref. 23, a mathematical model with Allee effect is proposed to simulate the impact of different dosages of RAI in the treatment of PTC patients. The model used was an ODE system, and the Allee threshold was associated with doses that may sufficiently eliminate tumors.

A mathematical model with two-step development is proposed in this article. The first was designed to examine the effectiveness of lenvatinib in the treatment of RAIR-DTC, in which the drug's action occurs by reducing the tumor carrying capacity, inhibiting tumor growth and altering the anti-immune capacity by activating NK cells. In the second step, we add the immunotherapeutic effect of pembrolizumab to the aforementioned dynamics, which, by preventing the binding b Drug developed and marketed by Merck & Co., Inc., under the label of Keytruda.

of PD-L1 to the PD-1 receptor, promotes the activation of antitumor T cells.

In what follows, Sections 2 and 3 introduce the model proposed, their parameters, and the analysis of the local and global stability of equilibria. Section 4 presents and discusses the results. Section 5 ends the article by articulating the major conclusions about the area of study.

Model with targeted therapy

We consider treatment with lenvatinib, with a concentration of the drug A, administered to patients with RAIR-DTC. Lenvatinib is administered orally in a daily dose of one to two capsules. [START_REF] Eisai | Highlights of prescribing information[END_REF] However, in our modeling we suppose that the concentration of the drug A is a continuous function in time. If the treatment cannot eradicate the tumor, it slows its growth to the maximum carrying capacity K 0 , this delay is represented by the variable S, which is a function of the concentration of the drug A. Assuming that N and L, respectively represent the number of cancer cells and NK cells, the model considered is the following ODE

                     dA dt = D -δA, dS dt = -c(A)S, dN dt = rN 1 - N K(S) -µAN -ρ(A)LN, dL dt = σ(A) -nL, (2.1)
where D, δ, r, µ, n are positive constants, the functions c, K : R + -→ R + are decreasing, and the functions ρ, σ : R + -→ R + are increasing. We set c 0 = c(0), K 0 = K(0), ρ 0 = ρ(0) and σ 0 = σ(0).

For the numerical simulations, we choose the functions

c(A) = c 0 e -cA , K(E) = K 0 e -KS , ρ(A) = ρ 0 + ρ 1 (1 -e -ρA )
and

σ(A) = σ 0 + σ 1 (1 -e -σA ),
where c 0 , c, K 0 , K, ρ 0 , ρ 1 , ρ, σ 0 , σ 1 and σ are positive constants, described in Table 1.

The functions ρ and σ simulate the increase in antitumor immune capacity during the treatment with lenvatinib, as proposed by Ref. 8. In the absence of treatment, D = 0, A → 0, the expressions used indicate a natural and limited response of the immune system, from the minimum values ρ 0 and σ 0 . The function σ incorporates the expected increase in the number of NK cells due to treatment, represented by the expression σ 1 (1 -e -σA ). The rate c simulates the velocity at which tumor cells reach the carrying capacity K 0 and, therefore, this rate should be decreasing as

The effect of Lenvatinib plus Pembrolizumab on thyroid cancer 5 a function of treatment with the antiangiogenic lenvatinib. The function K corresponds to the tumor growth towards carrying capacity, being directly affected by the antineoplastic action mentioned above through the control function S. Tumor cells can develop different mutations to resist to targeted therapy and continue to grow. This will be discussed later.

Equilibria and their local asymptotic stability

We consider the following parameter which is the equilibrium of the treatment

A d = D δ .
The trivial equilibrium, which always exists, is given by

E 0 := d, 0, 0, 1 n σ(d) .
A unique positive equilibrium is given by

E * := (d, 0, N * (d), L * (d)), with L * (d) = 1 n σ(d) and N * (d) = K 0 r r -µd - 1 n ρ(d)σ(d) .
A necessary and sufficient condition for the existence of the equilibrium E * is

f (d) := µd + 1 n ρ(d)σ(d) < r.
As the functions ρ and σ are positive and increasing, the function f is also positive and increasing with

f (0) = 1 n ρ 0 σ 0 and lim d→+∞ f (d) = +∞.
We have the following result.

• Suppose that 1 n ρ 0 σ 0 ≥ r. Then, f (d) ≥ r for all d ≥ 0. This means that only the trivial equilibrium E 0 = (d, 0, 0, σ(d)/n) exists.

• Suppose that 1 n ρ 0 σ 0 < r. Then, there exists a unique d 0 > 0 such that f (d 0 ) = r. In this case, if 0 ≤ d < d 0 , in addition to the trivial equilibrium E 0 , there exists a unique non-trivial equilibrium

E * := (d, 0, N * (d), L * (d)). If d ≥ d 0 , only the trivial equilibrium E 0 exists.
Let Ē := (d, 0, N (d), σ(d)/n) be any equilibrium of the system (2.1), with N = 0, for Ē = E 0 , and N = N * , for Ē = E * . The graphic of d → N (d) seen as a function of d is shown in Figure 1. We can see that N decreases and suddenly drops to 0 when d approaches d 0 . The four eigenvalues of the linearized system about the equilibrium point Ē are given by

-δ, -c(d), r 1 - 2 N (d) K 0 -f (d) and -n.
The local asymptotic stability of Ē depends only on the sign of the real eigenvalue 2.

λ(d) = r 1 - 2 N (d) K 0 -f (d).
• If 1 n ρ 0 σ 0 > r, then N (d) = 0 and λ := λ 0 (d) = r -f (d) < 0 for all d ≥ 0. In this situation, the trivial equilibrium E 0 is always locally asymptotically stable.

• Suppose that 1 n ρ 0 σ 0 < r and consider the unique

d 0 > 0 such that f (d 0 ) = r. If d > d 0 , the trivial equilibrium E 0 is the unique steady state with λ 0 (d) = r -f (d) < 0. Then, E 0 is locally asymptotically stable for d > d 0 . Now suppose that 0 ≤ d < d 0 . Then, the trivial equilibrium E 0 is unstable, because λ 0 (d) = r -f (d) > 0.
On the other hand, the positive equilibrium E * exists and the eigenvalue λ(d) := λ * (d) becomes

λ * (d) = -f (d) + r -2 r -µd - 1 n ρ(d)σ(d) = -r + f (d) < 0.
We conclude that E * is locally asymptotically stable for 0 ≤ d < d 0 .

The effect of Lenvatinib plus Pembrolizumab on thyroid cancer 7

:= [0, d] × R + × R + × [0, σ(d)
/n] is positively invariant under system (2.1). We make the change of variables

a = d -A and l = 1 n σ(d) -L.
Then, the system (2.1) becomes

                     da dt = -δa, dS dt = -c(d -a)S, dN dt = rN 1 - N K(S) -µN -ρ(d -a)(σ(d)/n -l)N, dl dt = σ(d) -σ(d -a) -nl.
(2.16)

As σ is an increasing function and a ∈ [0, d], we have

0 ≤ σ(d) -σ(d -a) = σ 1 e -σd [e σa -1] ≤ σ 1 σa.
For (a, S, N, l) ∈ Ω, we define the function

V 1 (a, S, N, l) = 1 2 (a 2 + S 2 ). The derivative of V 1 along the solution (a, S, N, l) ∈ Ω satisfies V1 (a, S, N, l) = -δa 2 -c(d -a)S 2 .
We proved that V 1 is a Lyapunov function on Ω. Let us define the set

Γ 1 := {(a, S, N, l) ∈ Ω : V1 (a, S, N, l) = 0}.
Using LaSalle's invariance theorem, we confirm that the set Γ 1 is globally attractive in Ω. Furthermore, for (a, S, N, l) ∈ Γ 1 , we have a = S = 0 and (N, l) satisfies the system

     dN dt = rN 1 - N K 0 -µdN -ρ(d)(σ(d)/n -l)N, dl dt = -nl.
(2.21)

Using the Lyapunov function V 2 (a, S, N, l) = 1 2 l 2 and LaSalle's invariance theorem on the set Γ 1 , we get the following globally attractive set of system (2.21)

Γ 2 = {(a, S, N, l) ∈ Γ 1 : V2 (a, S, N, l) = 0}.
For (a, S, N, l) ∈ Γ 2 , we have a = S = l = 0 and N satisfies the equation

dN dt = rN 1 - N K 0 -µdN - 1 n ρ(d)σ(d)N.
In other words

dN dt = r -f (d) - rN K 0 N. (2.24)
To conclude, we consider the two cases:

• If f (d) ≥ r, the trivial steady-state 0 of equation (2.24) is the only equilibrium and it is globally attractive for the solutions of (2.24). Then, the trivial equilibrium of system (2.1) is globally attractive on the set Ω.

• If f (d) < r, the positive equilibrium N * (d) = K 0 r (r -f (d))
of the logistic equation (2.24) is globally attractive. We conclude that the non-trivial equilibrium E * := (d, 0, N * (d), L * (d)) of system (2.1) is globally attractive.

In fact, one can prove the global attractivity of the two equilibria on the set (R + ) [START_REF] Nikiforov | Thyroid carcinoma: molecular pathways and therapeutic targets[END_REF] . The proof can be easily adapted from the above proof.

The antineoplastic action of lenvatinib as a retardant for the growth of the carrying capacity, expressed by the function c, is shown in Figure 2. The drug blocks factors as VEGF and its receptors, and causes the apoptosis of endothelial cells, among others. [START_REF] Rodia | Embolization of iliac metastasis during lenvatinib treatment in patient with advanced Hurthle cell thyroid carcinoma[END_REF] As seen, the growth of the tumor carrying capacity is not interrupted despite of treatment, but it evolves slowly over time, indicating a limited but important therapeutic outcome. The assigned values and units of each parameter are described in Table 2, which also presents the references used. A value was assigned to the parameter r based on the tumor doubling time observed in patients with metastatic PTC not responding to treatment with RAI 131 I, estimated in Ref. 26. For the parameters µ and ρ 1 , which represent important effects of the treatment, different scenarios with degrees of antitumor activity and drug efficacy are considered. In addition to the values of 1.4 ×10 -3 and 1 × 10 -10 , the values 1.5 × 10 -3 and 5 × 10 -10 are also used in the simulations. The parameter ρ 0 represents the immune response against the tumor in the absence of treatment. The parameter ρ, meanwhile, is associated with the increase in NK cells infiltrating the tumor, and thus an expansion of the antitumor immune response. The parameter σ 1 represents the constant activation rate of new NK cells due to the indirect action of lenvatinib. [START_REF] Zhang | Lenvatinib promotes antitumor immunity by enhancing the tumor infiltration and activation of NK cells[END_REF] The total number of cells activated in this manner represents an increase of approximately 1 × 10 7 over the expected number of cells responding naturally to the tumor, σ 0 . Considering that the rate σ 1 is positively affected by treatment with targeted therapy and that the concentration of lenvatinib stabilizes at 40.4 mg over time, [START_REF] Eisai | Highlights of prescribing information[END_REF] a value of 10 -1 is chosen for the parameter σ so that the equilibrium value for the population L is approximately 4 × 10 7 cells. In the next section, we include in the model (2.1) a new variable, P , which corresponds to the concentration of the immunotherapeutic agent pembrolizumab. The introduction of this antibody aims to promote the activation of antitumor T lymphocytes by preventing the binding of the PD-L1 protein to its PD-1 receptor.

Model with targeted therapy and immunotherapy

We consider in this section, a combined treatment between targeted therapy and immunotherapy in patients with RAIR-DTC, with the new variable P denoting the concentration of Pembrolizumab. We have the following model

                             dA dt = D -δA, dS dt = -c(A)S, dN dt = rN 1 - N K(S) -µAN -ρ(A)LN, dL dt = γ(A, P ) -nL, dP dt = Q -εP, (3.1)
where the new parameters Q and ε are positive, and instead of the function σ considered in the model ( 2), we introduce a new function γ : R + × R + -→ R + , which is assumed to be an increasing function of each component state. We set γ 0 = γ(0, 0) and for the numerical simulations, we take

γ(A, P ) = σ(A) + σ 2 (1 -e -σP ) = σ 0 + σ 1 (1 -e -σA ) + σ 2 (1 -e -σP ),
where the new parameters σ 2 and σ are positive. The new parameters are described in Table 3. Coefficient of activation of T cells by pembrolizumab

Local and global asymptotic stability of equilibria

As for the previous model, we denote

q = Q ε .
q is considered as a parameter for the second treatment. Similarly to the previously section, we obtain a trivial equilibrium E 0 := (d, 0, 0, γ(d, q)/n, q) which always exists. A unique non-trivial equilibrium is given by E * * = (d, 0, N * * (d, q), L * * (d, q), q), where

L * * (d, q) = 1 n γ(d, q) and N * * (d, q) = K 0 r r -µd - 1 n ρ(d)γ(d, q) .
A necessary and sufficient condition for the existence of E * * is given by

F (d, q) := µd + 1 n ρ(d)γ(d, q) < r.
For any fixed d ≥ 0, the function q → F (d, q) is increasing and satisfies

F (d, 0) = µd + 1 n ρ(d)σ(d) and F (d, +∞) = µd + 1 n ρ(d)(σ(d) + σ 2 ).
As a consequence, we obtain for a fixed d ≥ 0, the following properties.

• Suppose F (d, 0) > r. Then, for all q ≥ 0, F (d, q) > r and the trivial steady state E 0 is the only equilibrium. • Suppose F (d, +∞) < r. Then, for all q ≥ 0, F (d, q) < r and there exists a unique non-trivial equilibrium E * * .

• Suppose F (d, 0) < r < F (d, +∞). There exists a unique q 0 (d) > 0 such that F (d, q 0 (d)) = r. In this case, if 0 ≤ q < q 0 (d), there exists a unique non-trivial equilibrium E * * . If q ≥ q 0 (d), only the trivial equilibrium E 0 exists. Furthermore, d → q 0 (d) is a decreasing function.

As for the previous section, let Ê(d, q) = (d, 0, N (d, q), γ(d, q)/n, q) be any equilibrium of the system (3.1). From the linearized system, the five eigenvalues obtained are given by

-δ, -c(d), r 1 - 2 N (d, q) K 0 -F (d, q), -n and -ε.
Accordingly, the local stability of Ê(d, q) depends only of the eigenvalue

λ(d, q) = r 1 - 2 N (d, q) K 0 -F (d, q).
Define the set of existence of the non-trivial equilibrium E * * (d, q) R := {(d, q) : F (d, q) < r}.

It can be seen in Figure 3 that increasing the parameter d, which corresponds to treatment with lenvatinib, does not eradicate the tumor. Because, it is not possible to leave the region R of tumor stability. In particular, for a tumor that is not aggressive (r small), treatment with lenvatinib does not eliminate the tumor. On the other hand, the combination of the two treatments lenvatinib (d) and pembrolizumab (q) seems to be more effective, immunotherapy being the main responsible for the antitumor effect. Suppose that F (d, q) > r. Then, the only equilibrium is Ê(d, q) = E 0 , with N (d, q) = 0. This implies that λ(d, q) = r -F (d, q) < 0.

Hence, the trivial steady state E 0 is locally asymptotically stable. Now consider the case (d, q) ∈ R. Then, E 0 is unstable. On the other hand, the eigenvalue λ(d, q) associated to E * * (d, q) is given by

λ(d, q) = r -2 r -µd - 1 n ρ(d)γ(d, q) -F (d, q) = -r + 2F (d, q) -F (d, q) < 0.
We conclude that E * * (d, q) is locally asymptotically stable for (d, q) ∈ R.

As for system (2.1), using Lyapunov function and LaSalle's invariance theorem, one can prove the following result.

• If (d, q) / ∈ R, then the trivial steady state E 0 of system (3.1) is globally attractive.

• If (d, q) ∈ R, then the non-trivial equilibrium E * * (d, q) of system (3.1) is globally attractive.

For the numerical simulations, the values of the parameters of the model (3.1) are presented in Table 4. The parameter ε is calculated from the estimated half-life in Ref. 29 using the expression ε = ln [START_REF] Tirrò | Molecular alterations in thyroid cancer: from bench to clinical practice[END_REF] t 1/2 . The parameters σ and ¯σ are assumed to increase the number of NK cells and T lymphocytes to a total value of approximately 7 × 10 5 cells. This means that the total number of immune cells activated by the combination of drugs is more than three times that obtained without treatment. The ODE systems were solved using the fourth-order Runge-Kutta method, with a time step h = 0.01. In terms of modeling, only the situation with an attractive non-trivial equilibrium is confirmed by clinical data. The combination of treatments used does not usually eradicate the tumor. However, these treatments significantly slow the progression of the disease. 6, 30

Results and Discussion

The patients considered are RAIR-DTC, treated with lenvatinib A, or by combining lenvatinib A and pembrolizumab P . The effects of these therapeutic approaches include blocking the development of angiogenesis, represented by the function S and the carrying capacity K, the reduction in tumor growth rate, expressed as µAN , and the increase in the antitumor activity of the immune system with the activation of NK cells and T lymphocytes, as well as their infiltration into tumors, designated by the functions ρ, σ or γ.

The serum concentration limit for lenvatinib over time is approximately D/δ = 40.4 mg. This drug is administered orally in a daily dose of one to two capsules. [START_REF] Eisai | Highlights of prescribing information[END_REF] Pembrolizumab is an anticancer drug that is administered every 2-3 weeks with an elimination half-life of approximately 27.3 days. [START_REF] Du Rusquec | Clinical utility of pembrolizumab in the management of advanced solid tumors: an evidencebased review on the emerging new data[END_REF][START_REF] Mehnert | Safety and antitumor activity of the anti-PD-1 antibody pembrolizumab in patients with advanced, PD-L1-positive papillary or follicular thyroid cancer[END_REF][START_REF] Ahamadi | Model-based characterization of the pharmacokinetics of pembrolizumab: a humanized anti-PD-1 monoclonal antibody in advanced solid tumors[END_REF] In this study, the administration of pembrolizumab is considered daily with a dose of 15 mg. This is an approximation of a treatment given every 2 weeks at a dose of 2 mg/kg, as indicated in Ref. 29.

The data reported in Ref. 12 represent a study of 22 patients with advanced thyroid cancer treated with pembrolizumab for 24 months or until whether tumor progression or intolerable toxicity is confirmed, then in the numerical simulations the treatment period considered is 1000 days.

The natural influx of NK immune cells or T lymphocytes to the site of interaction with the tumor is considered in the numerical simulations with 3×10 5 cells. The use of lenvatinib generally activates NK cells, and the cellular increase corresponding to the use of this drug is shown in Figure 4, curve L(A). The combination of this treatment with pembrolizumab results in blockade of PD-L1 and PD-1 and therefore reduced the inhibition of the immune system by the tumor. [START_REF] Mehnert | Safety and antitumor activity of the anti-PD-1 antibody pembrolizumab in patients with advanced, PD-L1-positive papillary or follicular thyroid cancer[END_REF] In the curves representing the increase in immune cells, L(A) and L(A, P ), the initial number of immune cells, for t = 0, is assumed to be L(0) = 10 6 cells. Without any treatment, this cell population increases until a maximum of 3 × 10 7 cells. However, with the introduction of lenvatinib the maximum becomes approximately 4×10 7 cells. This slight increase is due to the fact that the expansion of the immune system observed in the experiments is a side effect not mentioned among the actions to be expected when using lenvatinib. The curve L(A, P ) represents the growth in the number of antitumor immune cells in terms of NK cells and T lymphocytes. With the addition of pembrolizumab, the maximum value becomes 10 × 10 7 cells, more than three times the number of immune cells in patients without treatment, curve L. The immunotherapeutic agent should cause a decrease in the tumor's ability to escape the immune response. However, in clinical practice, the approach has obtained only limited therapeutic effects. In Ref. 12, for example, among patients who showed partial response to treatment, the duration of the response ranged from 8 to 20 months. 2 and4. Initial conditions are S(0) = 1, A(0) = 24 mg, N (0) = 1.12 × 10 9 cells, L(0) = 10 6 cells and P (0) = 15 mg.

In general, the progression of iodine-refractory thyroid tumors under the treatments considered occurs in two ways: partial response, with a decrease in the number of cancer cells observed over a short period; and stable disease, in which tumor growth, despite an initial decrease, proceeds very slowly. [START_REF] Pitoia | Response to sorafenib treatment in advanced metastatic thyroid cancer[END_REF] Figure 5 illustrates the evolution of the tumor population with the two types of treatment: using only lenvatinib and using this drug in combination with pembrolizumab. Scenarios with different degrees of tumor elimination were obtained. In the simulations, although the initial decrease in the number of malignant cells is similar in both situations (t < 50 days), tumor evolution is different after the introduction of pembrolizumab, exhibiting greater capacity for removal of the tumor. Treatment with lenvatinib alone resulted in a partial tumor decrease within the first year, but the tumor started to grow again after this time. However, when considering the combination of treatments, the elimination of tumor cells is greater than that obtained with a single treatment, but the main result is that the tumor remains stable over time. 2 and4. Initial conditions are S(0) = 1, A(0) = 24 mg, N (0) = 1.12 × 10 9 cells, L(0) = 10 6 cells and P (0) = 15 mg.

The effect of the introduction of pembrolizumab in the treatment has a significant action on the tumor. In Figure 5, we consider some scenarios, by taking some values of the parameters µ and ρ 1 , of the action of the treatments on the behavior of the tumor. The parameter ρ 1 could be linked to the effectiveness of the two treatments, while µ refers only to the action of lenvatinib in the elimination of tumor cells. Assuming a change in the rate of tumor cell elimination due to the im-mune action, ρ 1 = 5 × 10 -10 , and without considering the immunotherapeutic drug pembrolizumab, Figure 6 shows a significant decrease in the number of malignant cells, but not enough to be similar to the combination of the two treatments. In the scenario assuming a greater efficacy of lenvatinib, µ = 1.5×10 -3 , there is a decrease of N , however it is less than that observed with the use of the two therapies. 2 and4, with the exception of the parameters µ and ρ 1 to which we give the values 1.5 × 10 -3 and 4 × 10 -10 , respectively. The initial conditions are S(0) = 1, A(0) = 24 mg, N (0) = 1.12 × 10 9 cells, L(0) = 10 6 cells, and P (0) = 15 mg.

Similar to lenvatinib, sorafenib is a multiple kinase inhibitor (MKI) used for different types of malignant tumors. [START_REF] Roskoski | Properties of FDA-approved small molecule protein kinase inhibitors[END_REF] The observational study by Ref. 31 examined the efficacy of sorafenib in the treatment of metastatic thyroid carcinoma RAIR. In this study, 8 patients received 400 mg of the drug per day until disease progression or intolerable toxicity was observed. One patient showed partial response, with 35% tumor regression after 6 months of treatment, while 5 others showed stable disease. In the review conducted by Ref. 33, the authors found that the development of resistance after 1 or 2 years of treatment may occur in patients with a partial response or with stable disease. They then recommend using different MKIs. Various treatment options for RAIR-DTC are discussed in Ref. 15, which highlights the poor performance of some conventional treatments, including administration of chemotherapy via doxorubicin, a treatment recommended only when MKIs are not indicated. In discussing the use of lenvatinib to treat anaplastic thyroid cancer, Ref. 34 identified resistance mechanisms that could be activated by tumors to survive treatment. RNA-seq analyses from this study revealed that prolonged treatment with lenvatinib caused changes in several cellular and oncogenic pathways. In another study of 94 patients regarding the efficacy and safety of lenvatinib as a treatment for RAIR-DTC, Ref. 35 observed partial responses and stable disease in only 36% and 41% of patients, respectively.

The antiangiogenic effect of lenvatinib, represented in our study by tumor carrying capacity, K, aims to mimic the action of the drug by blocking certain factors, for example VEGF or its receptor, and causes apoptosis of endothelial cells. [START_REF] Rodia | Embolization of iliac metastasis during lenvatinib treatment in patient with advanced Hurthle cell thyroid carcinoma[END_REF] However, since the treatment has a limited effect, we propose that the tumor develops new molecular changes which re-expand its carrying capacity and reactivate its tumor growth process. [START_REF] Khan | Targeting XPO1 and PAK4 in 8505C anaplastic thyroid cancer cells: putative implications for overcoming lenvatinib therapy resistance[END_REF] The action of lenvatinib as an inhibitor of fibroblast growth factor (FGF) and platelet-derived growth factor (PDGF), involved in the processes of growth and cell division, is also represented by the term µAN .

For the numerical simulations, several values of the parameters representing the efficacy of lenvatinib against tumor growth µ and the rate of elimination of malignant cells via the immune action ρ 1 , also enhanced by the second drug, were tested. We changed only these two parameters because they represent the direct action of drugs. Therefore, by considering doses similar to those taken in clinical studies, our model produces results that should expand current knowledge about targeted therapy and immunotherapies applied to RAIR-DTC.

Conclusions

A two-step development ODE model have been proposed to test therapies applied to patients with RAIR-DTC treated with lenvatinib alone or in combination with pembrolizumab. Numerical simulations have shown the effectiveness of the treatments in inducing a partial response or a stable disease. This response is more effective in the case of a combination of the two therapies. The parameters related to the action of the drugs were analysed and scenarios with different responses to the treatment were obtained. Our study provides new ideas for treatment discussions for patients with RAIR-DTC and highlighted the importance of combination therapies, including targeted therapies and immunotherapies.

Fig. 1 .

 1 Fig. 1. Evolution of N as a function of d, the other parameters are chosen from Table2.

Fig. 2 .

 2 Fig. 2. Progression of tumor carrying capacity in RAIR-DTC patients treated or not with lenvatinib. Parameters used as shown in Table2and initial conditions S(0) = 1, K 0 = 10 10 cells and D = 0, 6, 12 and 24 mg.

Fig. 3 .

 3 Fig.3. Contour levels F (d, q) = r for some values of r in the plane d × q. Parameters used as shown in Table2with the exception of µ = 4 × 10 -10 and ρ 0 = 1 × 10 -9 .

Fig. 4 .

 4 Fig. 4. Evolution of cell populations of NK and T lymphocytes without therapy, L, and during treatment with lenvatinib, L(A), and the combination of lenvatinib with pembrolizumab, L(A, P ). Parameters are taken from Tables2 and 4. Initial conditions are S(0) = 1, A(0) = 24 mg, N (0) = 1.12 × 10 9 cells, L(0) = 10 6 cells and P (0) = 15 mg.

  ], µ = 0.0014, ρ 1 = 10 -10 N [ A, P ], µ = 0.0014, ρ 1 = 10 -10

Fig. 5 .

 5 Fig. 5. Evolution of the population of malignant cells during treatments with lenvatinib alone, N [A], and lenvatinib combined with pembrolizumab, N [A, P ]. Parameters are taken from Tables2 and 4. Initial conditions are S(0) = 1, A(0) = 24 mg, N (0) = 1.12 × 10 9 cells, L(0) = 10 6 cells and P (0) = 15 mg.

10 Fig. 6 .

 106 Fig. 6. Evolution of the malignant cell population during treatment with lenvatinib only, N [A] and in combination with pembrolizumab, N [A, P ]. The parameters are taken from Tables2 and 4, with the exception of the parameters µ and ρ 1 to which we give the values 1.5 × 10 -3 and 4 × 10 -10 , respectively. The initial conditions are S(0) = 1, A(0) = 24 mg, N (0) = 1.12 × 10 9 cells, L(0) = 10 6 cells, and P (0) = 15 mg.

Table 1 .

 1 Parameters and their meanings.

		Description
	D	Recommended daily dose of lenvatinib
	δ	Elimination rate of lenvatinib drug
	c 0	Minimum rate at which tumor cells reach the carrying capacity
	c	Coefficient of the antineoplastic effect of lenvatinib
	r	Proliferation rate of tumor cells
	K 0 Maximum tumor carrying capacity
	K	Coefficient of growth of tumor carrying capacity
	µ	Rate of drug efficacy in eliminating tumor cells
	ρ 0	Elimination rate of tumor cells due to immune system action
	ρ 1	Elimination rate of tumor cells due to immune action with lenvatinib
	ρ	Coefficient of immune action due to lenvatinib
	σ

0 Coefficient of the natural influx of immune cells to the tumor site σ 1 Limit rate of activation of immune cells due to lenvatinib σ Rate of activation of NK cells by lenvatinib n Natural mortality rate of immune cells 2.2. Global asymptotic stability of the equilibria It is clear that the set Ω

Table 2 .

 2 Parameters with values, units and references.

	March 8, 2022 14:51	"Published paper"
	10 J.G. Da Silva, I.C. Da Silva, M. Adimy & P.F. Mancera
		Parameter Value	Unity	Reference
		D	24	mg × day -1	Ref. 6
		δ	0.5939	day -1	Ref. 24
		c 0	10 -2	day -1	assumed
		c	5 × 10 -2	mg -1	assumed
		r	0.0707	day -1	Ref. 26
		K 0	10 10	cells	Ref. 27
		K	1	-	assumed
		µ	1.4 × 10 -3 (mg × day) -1	assumed
		ρ 0	1 × 10 -10	(cells × day) -1 assumed
		ρ 1	1 × 10 -10	(cells × day) -1 assumed
		ρ	10 -2	-	assumed
		σ 0	3 × 10 5	cells × day -1	Ref. 28
		σ 1	1 × 10 5	cells × day -1	assumed
		σ	10 -1	-	assumed
		n	10 -2	day -1	assumed

Table 3 .

 3 Meaning of the new parameters. Limit rate of activation of T cells due to the action of pembrolizumab σ

		Description
	Q Recommended daily therapeutic dose of pembrolizumab
	ε	Elimination rate of pembrolizumab
	σ 2

Table 4 .

 4 Parameters with values, units and references.

	Parameter Value	Unity	Reference
	Q	15	mg × day -1	Ref. 29
	ε	0.0254	day -1	Ref. 12
	σ 2	6 × 10 5 cells × day -1 assumed
	σ	10 -1	-	assumed

Acknowledgments

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brasil (CAPES) -Finance Code 001.


Modeling time resolved photoluminescence

The model we used to interpret the decays is the 1d-Drift-Diffusion model. This model is the

following. (Ref S1)
We consider a slab of intrinsic semiconductor of thickness 𝐿 (the thickness is changed in the model for each etching time). The laser pulse will generate excited charges in the device. We consider that the electron and hole density to be equal at all position in the thickness and time.

We model the time and space (in thickness) dependent photo generated carrier density Δ𝑛(𝑧, 𝑡)

with the 1d-Drift-Diffusion equation:

𝜕Δ𝑛 𝜕𝑡 = 𝐷 𝜕 2 Δ𝑛 𝜕𝑧 2 -𝑘 1 Δ𝑛 -𝑘 2 Δ𝑛 2 + 𝑔(𝑧, 𝑡) (S1)
With the effective diffusion coefficient 𝐷, the SRH recombination constant 𝑘 1 and the radiative external recombination coefficient 𝑘 2 . The laser pulse is modeled via the time dependent generation rate 𝑔 via:

𝑔(𝑧, 𝑡) = [𝑛 𝛾 𝛼𝑒 -𝛼𝑧 ] × [ 1 √2𝜋𝜎 exp ( -𝑡 2 2𝜎 2 )]
That is Beer-Lambert's absorption law with 𝛼 the absorption coefficient at the laser wavelength 𝛼(532𝑛𝑚) = 5e4 𝑐𝑚 -1 and 𝑛 𝛾 the fluence in ph.cm -2 pulse -1 which is varied during the experiment. The right part of 𝑔 is a Gaussian temporal profile of duration 𝜎 = 5𝑝𝑠 in the simulation.

The boundary conditions to solve Eq S1 are the following:

{ 𝜕Δ𝑛 𝜕𝑧 | 𝑧=0 = 𝑆 𝑡𝑜𝑝 Δ𝑛(𝑧 = 0, 𝑡) 𝜕Δ𝑛 𝜕𝑧 | 𝑧=𝐿 = -𝑆 𝑏𝑜𝑡 Δ𝑛(𝑧 = 𝐿, 𝑡)
They represent non-radiative recombination at the interfaces.

The implementation of the model is made via a Matlab code using the pdepe function to solve the partial differential equation. A fitting procedure was coded and has the following properties:

1. A unique model is fitted for all the selected fluences. That means that only one value for the physical parameters is fittedno fluence dependent values are fitted.

2. The uncertainties on the fitted parameters are estimated using the nlparci method. S3. Model parameters values.

The model parameters are presented in Table S1. We chose to estimate the bulk recombination coefficient 𝑘 1 to 10 6 𝑠 -1 as this value is close to what we found previously for our baseline perovskite.

To analyze the impact of the etching on the recombination, we fitted for each set of curves (ad of Fig. S2) the drift diffusion model with only one varying parameter: the top surface recombination velocity 𝑆 𝑡𝑜𝑝 . The results are given in Fig S2 (e). We observe that the top surface recombination velocity is found to decrease with the decreasing thickness. This suggests that the etching did not compromise the quality of the surface. With the fitted parameters, bulk nonradiative recombination account for 22% (400nm-thick device) to 27% (150nm-thick device) of all recombination at 1 sun equivalent illumination, while top-surface recombination account for 77% (400nm-thick device) to 71% (150nm-thick device). The remaining are radiative recombination.

Impact of the value of the parameter k1 on fitting results

Modeling the effect of thickness on PL peak

In order to decorrelate the effect of thickness and potential surface modification on steady-state hyperspectral PL acquisitions, we performed drift-diffusion simulation with Atlas solver. We modelled the effect of thickness decrease keeping material properties constant in order to conclude on its impact on PL. We observed in Figure S3 that when the thickness is decreased, the PL peak decreases which is in agreement with the observed analysis. When the thickness decreases, the surface recombination exhibits a more significant effect. The surface/volume ratio increases which leads to a decrease in bulk recombination versus surface recombination.

Complementary experimental analysis detailed in the main section of the manuscript aids in deducing whether the thickness is the only factor impacting the PL peak intensity decrease or if the surface recombination is changed at the same time.

Figure S1 .

 S1 Figure S1. XPS high-resolution spectra of Br 3d, C 1s, N 1s, I 3d, Pb 4f and Cs 3d for the 4 studied samples. The C -N contribution at 288.3 ± 0.2 eV corresponds to the FA + environment and the peak at 286.4 ± 0.2 eV to the MA + one. The obtained spectra for the 4 samples are

Figure S2 .Figure S3 :

 S2S3 Figure S2. XPS high-resolution spectra of Pb 4f, N 1s, C 1s, I 3d and Br 3d, before sputtering

Figure S5 .

 S5 Figure S5. Influence of the perovskite thin film thickness on the maximum PL intensity. The incident illumination power was varied from 1 W/cm 2 to 1e5 W/cm 2 .

Figure S6 .

 S6 Figure S6. PL average spectra and corresponding fitting curve for reference sample, crater 2

Figure S8 :

 S8 Figure S8: Hyperspectral acquisition in different location of the crater. (left) Corresponding

Table S1 .

 S1 Atomic

	Crater number Estimated remaining	Average Roughness
			thickness [nm]	(RMS) [nm]
	Surface before		400	14 ± 2
	sputtering
	2		250	17 ± 2
	3		200	17 ± 3
	4		150	18 ± 3
	5		100	16 ± 4
	6	Interface perovskite	20 ± 5
			/TiO2
	Atomic %	Fresh Perovskite After 1h analysis
	I3d5/2		30.9	30.7
	Pb4f7/2		9.5	10.0
	N1s		18.7	18.7
	Cs3d5/2	0.6	0.5
	C1s		33.9	33.8
	O1s		0.6	0.8
	Br3d		5.8	5.7

percentage of I3d5/2, Pb4f7/2, N1s, Cs3d5/2, C1s, O1s and Br3d for fresh Perovskite and after 1h of analysis under X-ray beam.

Table S2 .

 S2 Values

of the average roughness measured via Gwyddion software for different estimated thicknesses in the perovskite layer.

Table

	Bulk defect SRH	𝑘 1	1𝑒6 𝑠 -1	Estimation from
	coefficient			previous knowledge
	Radiative	𝑘 2	3 × 10 -10 𝑐𝑚 3 𝑠 -1	Fitted from the full
	recombination			cell reference
	coefficient
	Diffusion coefficient	𝐷	1 × 10 -2 𝑐𝑚 2 𝑠 -1	From a fit of the thin
				film reference.
	Bottom surface	𝑆 𝑏𝑜𝑡	0𝑐𝑚. -1	Hypothesis to
	recombination			compare the impact
	velocity			of the top surface
	Absorption	𝛼(532𝑛𝑚)	5e4 𝑐𝑚 -1	Estimated from
	coefficient at 532nm			absorption
				measurements
				They
	are an estimate of the uncertainties -which are difficult to determine for this non-linear
	fitting method. (Ref S2)
	Parameter	Symbol	Value	Comment
	Thickness	𝐿	Varied for each set of	From the expectation
			curves

Supporting Information

Indepth chemical and optoelectronic analysis of triple-cation Perovskite thin films by combining

XPS profiling and PL Imaging

Stefania Cacovich † (1) , Pia Dally † (2,3) , Guillaume Vidon (2) , Marie Legrand (2,4) , Stéphanie Gbegnon (2) , Jean Rousset (2,4) , Jean-Baptiste Puel (2,4) , Jean-François Guillemoles (1) , Philip Schulz (1) , Muriel Bouttemy (2,3) , Arnaud Etcheberry (3) To check the influence of the parameter k1 we reperform the fit with different values of k1.

The base value was set to k1=1e6 s-1. We test two other values: dividing by two or multiplying by two. We obtain the following fitting results: We see that some conclusions are true irrespective of the value of k1: the value of the surface recombination parameter decreases for thinner craters. This confirms that the Ar Etching is not causing major damage to the surface recombination in the TR-PL experiment.

When the value of k1 is reduced (orange bars), the surface recombination is increased which is coherent: the model needs to have a certain amount of non-radiative recombination that can either come from the surface or the bulk of the material. As we mentioned in the main text, this makes that the parameters k1 and Stop are highly anti correlated. In [START_REF] Hutter | Pitfalls and prospects of optical spectroscopy to characterize perovskite-transport layer interfaces[END_REF], Hutter et al explain how TR-PL is a powerful tool but cannot distinguish different sources of non-radiative recombination one from the other. That is why we choose to fix the value of k1 and to interpret the variations of non-radiative recombination in terms of surface recombination.