Johannes Wirtz
email: jwirtz@lirmm.fr

Stéphane Guindon
email: guindon@lirmm.fr

On the connections between the spatial Lambda-Fleming-Viot model and other processes for analysing geo-referenced genetic data

Keywords: Spatial Lambda-Fleming-Viot, Birth-Death processes, Duality, Efficient simulation

. The ΛV model provides a sound mathematical framework for describing the evolution of a population of related individuals along a spatial continuum. It alleviates the "pain in the torus" issue with Wright and Malécot's isolation by distance model and is sampling consistent, making it a tool of choice for statistical inference. Yet, little is known about the potential connections between the ΛV and other stochastic processes generating trees and the spatial coordinates along the corresponding lineages. This work focuses on a version of the ΛV whereby lineages move infinitely rapidly over infinitely small distances. Using simulations, we show that the induced ΛV tree-generating process is well approximated by a birth-death model. Our results also indicate that Brownian motions modelling the movements of lineages along birth-death trees do not generally provide a good approximation of the ΛV due to habitat boundaries effects that play an increasingly important role in the long run. Finally, we describe efficient algorithms for fast simulation of the backward and forward in time versions of the ΛV model.

Introduction

The integrated analysis of genetic and spatial data in the fields of phylogeography or spatial population genetics is central to our understanding of the forces driving the evolution of living organisms in space and time. Indeed, accommodating for the evolutionary relationships between individuals of the same population or between distantly related species when analysing their spatial distribution permits the reconstruction of ancestral migration and dispersal events. It then becomes possible to examine the links between these events and past environmental or ecological changes so as to decipher the mechanisms underlying key biological processes such as speciation or the impact of natural selection on spatial patterns of biodiversity. Combining the horizontal (spatial) and vertical (evolutionary) dimensions of geo-referenced genetic data is therefore paramount in order to elucidate the mechanisms and test hypotheses about the underlying data generating processes.

In population genetics, Wright's island model [START_REF] Wright | Evolution in Mendelian populations[END_REF] was the first of a series of "migration-matrix" models that aimed at describing the evolution of a population that is spatially structured in distinct demes (see Rousset [START_REF] Rousset | Inferences from spatial population genetics[END_REF] for a review). Despite their relative simplicity, the island model and its descendants, including most notably the stepping stone model [START_REF] Kimura | Stepping stone' model of population[END_REF], provided population geneticists with a rich set of tools to test important biological hypotheses such as panmixia or the existence of past and/or ongoing migrations between sub-populations.

The assumption of discrete demes is convenient mathematically. The ability to accommodate for populations that are spatially distributed along a continuum is a natural extension of the discrete assumption. That extension is expected to significantly expand the range of applications and, in numerous instances, enhance the relevance of spatial population genetics models [START_REF] Bradburd | Spatial population genetics: it's about time[END_REF]. Over the last eight decades, progresses in the development of these models turned out to be rather slow and faced serious difficulties in some cases. The isolation by distance model proposed by Sewall Wright [START_REF] Wright | Isolation by distance[END_REF] and Gustave Malécot [START_REF] Malécot | Mathematics of heredity[END_REF], for instance, was shown to suffer from pathological behaviour in the long run (the so-called "pain in the torus" described by Joseph Felsenstein,[START_REF] Felsenstein | A pain in the torus: some difficulties with models of isolation by distance[END_REF]), forcing population geneticists to rely on the discrete approximation aforementioned. The approaches proposed by Wilkins and Wakeley [START_REF] Wilkins | The coalescent in a continuous, finite, linear population[END_REF], Wilkins [START_REF] Wilkins | A separation-of-timescales approach to the coalescent in a continuous population[END_REF] addressed the "clumping" issue that hampered the isolation by distance model. Yet, as pointed by Alison Etheridge and colleagues, the models proposed here lacked sampling consistency, implying that the time to coalescence of lineages depended on the size of the sample considered [START_REF] Barton | A new model for evolution in a spatial continuum[END_REF], thereby limiting their application in practise.

While there are relevant approaches available that provide graphical summaries of populations distributed along a spatial continuum (see e.g., [START_REF] Novembre | Genes mirror geography within europe[END_REF][START_REF] Wang | A quantitative comparison of the similarity between genes and geography in worldwide human populations[END_REF][START_REF] Bradburd | A spatial framework for understanding population structure and admixture[END_REF][START_REF] Bradburd | Spatial population genetics: it's about time[END_REF]), sound mechanistic models that accommodate for continuous diffusion of individuals in their habitat along with genetic drift are scarce. In a pioneering work, Alison Etheridge, Nick Barton and Amandine Véber [START_REF] Barton | A new model for evolution in a spatial continuum[END_REF][START_REF] Barton | Modelling evolution in a spatial continuum[END_REF] introduced the spatial Lambda-Fleming-Viot model (noted as ΛV in the following) in an attempt to fill this gap. To the best of our knowledge, the ΛV is the sole mechanistic model that (1) accommodates for populations distributed along a spatial continuum, under a stationary regime (i.e., the population density does not change, on average, during the course of evolution) and (2) provides a coherent account of the forward in time evolution of a population along with a dual description of the backward in time evolutionary dynamics of a sample from that population and (3) is amenable to parameter inference using a Bayesian approach, applicable to small to moderate size data sets, e.g., see [START_REF] Guindon | Demographic inference under the coalescent in a spatial continuum[END_REF][START_REF] Joseph | Demographic inference under a spatially continuous coalescent model[END_REF]).

The properties of the ΛV model and some extensions are well characterised mathematically [START_REF] Véber | The spatial Lambda-Fleming-Viot process: An event-based construction and a lookdown representation[END_REF][START_REF] Biswas | The spatial lambdafleming-viot process with fluctuating selection[END_REF][START_REF] Louvet | Stochastic measure-valued models for populations expanding in a continuum[END_REF]]. Yet, relatively little is known about the relationships between ΛV and other popular population genetics models. Shedding light on potential connections between these models would help delineate conditions in which the ΛV may be well approximated by other processes, potentially leading to more efficient parameter estimation procedures. More importantly, establishing such bridges would help gain a better understanding of the biological relevance of the ΛV process.

In this study, we consider the non-trivial case where the rate of reproduc-tion and extinction (REX) events in the ΛV model is large and the radius of each event (i.e., the parent-to-offspring distance) is small. We first focus on the tree-generating process that derives from the ΛV model forward in time in these particular conditions and show that the distribution of trees deriving from the ΛV is well approximated by that obtained from a birth and death (BD) process. We then incorporate the spatial component in our analyses and show how the ΛV model compares to the birth and death model with spatial coordinates fluctuating along lineages according to a Brownian process, as introduced in [START_REF] Lemey | Phylogeography takes a relaxed random walk in continuous space and time[END_REF] and available in the popular software package BEAST [START_REF] Drummond | BEAST: Bayesian evolutionary analysis by sampling trees[END_REF][START_REF] Suchard | Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10[END_REF]. Results from simulations indicate that habitat border effects that come into play with the ΛV model but are ignored by the Brownian process, preclude the convergence of both models to the same process. Finally, we describe two algorithms for efficient simulation of the ΛV process forward and backward in time, which are at the core of some of the model comparisons performed here.

Notation and models

We first introduce some notation that will be used throughout the manuscript. Let n be the number of sampled lineages. τ denotes a ranked tree topology with n tips and t, the corresponding vector of 2n -1 node times, which are defined relative to the sampling time. Throughout this study, sampling of lineages takes place at a single point in time (i.e., we do not account for heterochronous data) taken to be equal to 0. ℓ is the vector of 2n -1 spatial coordinates at all nodes in the tree.

Individual-based ΛV model

We consider the forward-in-time version of this process here, taking place on a w×h rectangle, denoted as A in what follows. Individuals that constitute the population of interest are distributed uniformly at random with density ρ on that rectangle. Lineage reproduction and extinction (REX) events occur at rate ξ, the per unit space rate. When one such event takes place, (1) individuals die with probability υ exp(-d 2 /2θ 2 ), where d is the Euclidean distance between the corresponding individual position and the location of the center of the REX event, (2) offspring are generated according to a nonhomogeneous Poisson process with intensity ρυ exp(-d 2 /2θ 2 ) and (3) one parent for the newly generated offspring is chosen where a parent at distance d from the centre has probability proportional to exp(-d 2 /2θ 2 ) to be selected (individuals that die on that event may also be selected as parent).

A closed-form formula for the likelihood, i.e., the joint probability density of τ and t conditioned on υ, θ, ξ and ρ is not available. Yet, obtaining random draws from the corresponding distribution is relatively straightforward. In particular, in a manner similar to the Wright-Fisher model and Kingman's coalescent [START_REF] Kingman | On the genealogy of large populations[END_REF], the ΛV has a backward in time dual of the forward in time process that allows for rapid simulations of genealogies of a sample of n lineages (see [START_REF] Barton | A new model for evolution in a spatial continuum[END_REF] and section 2.5.2 for a description of an efficient backward in time algorithm for simulating a two-tip genealogy under the ΛV).

Birth and death process with Brownian diffusion

Beside the ΛV model, this study focuses on the homogeneous BD model with complete sampling. According to this process, a first lineage arises at time t or , the time of origin of the process. The rate at which any given lineage splits/dies is governed by the birth and death parameters λ and µ respectively, i.e., the process is homogeneous so that per-lineage birth and death rates are fixed throughout. Data collection takes place in the future compared to t or , at which point we condition the genealogy τ on having n live lineages, which are all included in our sample. Spatial coordinates evolve in a two dimensional space, with movements along the northings considered as independent from that along the eastings. In each dimension, the spatial position of a lineage fluctuates according to a Brownian process with diffusion parameter σ. Hence, the distribution of the position at the end of a branch of length t (in calendar time units) is Gaussian with mean given by the lineage position at the start of that branch and variance σ 2 t.

The very idea of using Brownian diffusion to model the evolution of locations along a genealogy was introduced in [START_REF] Lemmon | A likelihood framework for estimating phylogeographic history on a continuous landscape[END_REF][START_REF] Lemey | Phylogeography takes a relaxed random walk in continuous space and time[END_REF]. Although Lemey et al. [START_REF] Lemey | Phylogeography takes a relaxed random walk in continuous space and time[END_REF] focused on a "relaxed" version of this approach, whereby each branch in the phylogeny has its own spatial diffusion parameter, we focus here instead on the "strict" version of the model, with a single diffusion parameter applying to all edges of the tree. This model is noted BD 2 in the following for BD model combined with Brownian Diffusion.

Connecting the two models

Our study first focuses on the comparison between between p BD (τ, t|λ, µ, n, σ), the likelihood of the BD model, and the equivalent density for the ΛV model, p ΛV (τ, t|ξ, θ, n, υ, ρ), when focusing only on the tree-generating parts of both models. We then examine the link between p BD 2 (τ, t, ℓ|λ, µ,n, σ) and p ΛV (τ, t, ℓ|ξ, θ,n, υ, ρ), the full likelihoods, i.e., including both the tree and the spatial components, of ΛV and BD 2 . We consider the particular case where ξ → ∞ and θ → 0, i.e., REX events occur at a high rate and each of them has a very small radius. We assume here that ξθ 2 → c for some c ∈ R. Disregarding the spatial component, we refer to the "limit" model as ΛV * .

Birth and death rates under ΛV

We first focus on the rate µ ΛV (x) at which an individual at position x ∈ A dies in a REX event under ΛV. Events occur uniformly on A at rate ξ, and given an event location z ∈ A, the probability that an individual at position x dies due to the event is υ exp(-d 2 /2θ 2 ), where d = ∥x -z∥. So the overall rate at which an individual at x dies is obtained by integrating this probability over the habitat, multiplied by ξ|A|, where |A| is the area of the habitat. So we have:

µ ΛV (x) = ξ|A| A υ |A| exp(-∥x -z∥ 2 /2θ 2 )dz (1)
When θ → 0, ξ → ∞ and θ 2 ξ → c, the right-hand side can be written as

µ ΛV * (x) = 2πcυ, (2)
We observe that the rate hence obtained does not depend on the individual's position x. In particular, this rate does not depend on the distance to the edges of the habitat. Therefore, in the ΛV * all individuals have a unique "death rate" µ ΛV * = 2πcυ.

Similarly, individuals have a rate of "giving birth" analogous to the birth rate of a BD process. At any REX event, the cumulative intensity of births on A is

A ρυ exp(-∥y -z∥ 2 /2θ 2 )dz → 2πθ 2 ρυ (3)
where z denotes the event location. Therefore, in the ΛV * process, the number of individuals born in one event is Poisson with parameter 2πθ 2 ρυ and the probability that in one event k individuals are born is

p k := (2πθ 2 ρυ) k k! exp(-2πθ 2 ρυ) (4)
In [START_REF] Barton | Modelling evolution in a spatial continuum[END_REF], it is shown that the probability that an individual at location x is chosen as the parent by an event located at z is

1 2πθ 2 ρυ exp(-∥x -z∥ 2 /2θ 2 ) • (1 + O(ρ -1 )) (5)
We shall assume that ρ is large enough such that the order term in (5) becomes negligible. However, we note that when simulating we observed that even for values of ρ around 1 this seemed to be the case on average. Combining ( 5) and ( 4), we can calculate the rate at which an individual at x is chosen as the parent by a REX event and k individuals are being generated by that event as

λ (k) ΛV (x) = ξ|A| A p k • 1 2πθ 2 ρυ|A| exp(-∥x -z∥ 2 /2θ 2 )dz (6)
Now, letting ξ → ∞ tend to infinity and θ → 0 in the same way as before, we have lim λ

(1)

ΛV (x) = lim ξ ρυ p 1 = 2πc (7) and λ (k)
ΛV (x) → 0 for all k > 1. The limit thus eliminates the possibility of multiple offspring during one event, ensuring that an individual can give birth to at most one child at a time. The rate at which lineages split in the ΛV * is then

λ ΛV * (x) = 2πc (8)
In the standard BD model, birth and death events never happen at the very same point in time, let alone along the same lineage. In the ΛV model however, a REX event may involve the splitting of the parental lineage and the death of that same lineage. Given a REX event with centre z, the probability for a given lineage located at x to die or to give birth to new lineages is proportional to exp (-∥x -z∥ 2 /2θ 2 ). The probability of both events (birth and death) taking place is thus proportional to exp (-∥x -z∥ 2 /θ 2 ). Birth or death events alone therefore become infinitely more probable than simultaneous birth and death events when the radius tends to zero so that, in that respect, the ΛV behaves in a manner similar to the BD model.

We conclude that in the ΛV, for diminishing values of θ and increasing ξ, the number of offspring lineages is stochastically similar to the number of lineages in a birth and death process with λ = 2πc, µ = 2πcυ, where c = θ 2 ξ is constant. Since 0 < υ ≤ 1, we always have λ ≥ µ, so the process is supercritical, except in the case where υ = 1. We shall make use of the latter assumption throughout this manuscript.

We confirmed these observations by simulating the ΛV forward in time. At the beginning of each simulation run, we randomly selected one individual within the population. Simulations stopped whenever this individual was the target of an event, and the time at which that event took place was recorded. Three types of events can be observed: 1) The death of the individual; 2) the individual giving birth to one or more offspring individuals; and 3) death and birth of that individual at the same time. The values for θ and ξ were chosen in such a way that c = θ 2 ξ = 1 throughout our simulations. Also we set υ = 1, ρ = 20, w = 10 and h = 10.

As θ decreases (and ξ increases), we observe that events of the third type (simultaneous birth and death of the same lineage) become increasingly rare, while events of type one (death) and two (birth) occur at about the same frequency. For example, for θ = ξ = 1 the frequency of type three events observed is close to 0.08, rises to 0.15 for θ = 0.25, ξ = 16, then drops to 0.01 for θ = 0.05 and to effectively zero for smaller θ. For type one, the frequencies are 0.74, 0.63, 0.54 and 0.50, whereas the frequencies of events of type three are 0.18, 0.22, 0.45 and finally about 0.50 as well.

From the times recorded at which these events take place, we reconstruct the probability density of the time to an event of the respective type. These densities are represented in Figure 1 for various values of θ, while the black curves represents the densities derived for death events (right) and birth events (left), respectively in a BD process with parameter λ = µ = 2π, which both conform to an exponential density with parameter λ = 2π. For decreasing θ, we observe a trend of the densities in the ΛV to approach those in the BD.

On the offspring number distribution in the two processes

In this section, we examine the distribution of the number of descendant lineages resulting from each individual in the initial population with respect to time. More specifically, we monitor the number of live descendants of every individual in the starting population. The number of descendants of each of these ancestors at some time t defines its family size. In the BD process, since all individuals are independent, the evolution of the size of a family behaves the same way as the size of a population (i.e., the number of surviving lineages) that started with a single individual.

When focusing on the fate of a single ancestor, both the BD and the ΛV processes have one absorbing state: Whenever a family size reaches zero, the processes stay in that state (the family has become "extinct"). The BD processes that correspond to ΛV processes are either critical (λ = µ, which is the case we consider here) or supercritical (λ > µ). In the critical BD, when starting with one individual at time 0, the process will eventually reach 0 with probability one, i.e., any family will become extinct after a sufficient amount of time (although the expected time to that event is infinite). Starting from one individual at time 0 and conditioning on non-extinction, the probability p * 1m (t) of observing a certain family size m > 0 at time t for the critical process is given by

p * 1m (t) = (λt) m-1 (1 + λt) m (9)
as stated in [START_REF] Tavaré | The linear birth-death process: an inferential retrospective[END_REF] (Equation 4). It is noteworthy that p * 1m (t) converges in distribution if and only if the process is subcritical (λ < µ) [START_REF] Karlin | Chapter 9 -stationary processes[END_REF][START_REF] Cavender | Quasi-stationary distributions of birth-and-death processes[END_REF], which the ΛV is unable to emulate.

Since birth and death rates in the ΛV* correspond to a critical BD, the family size of any individual from the initial population evolving under ΛV is expected to drop to 0 after some (potentially infinite) time. On the other hand, while the death rates are constant and the same for all individuals in the ΛV, the birth rate of one individual may be affected by the number of individuals close by; for example, if a neighbourhood C ⊆ A is momentarily sparsely populated, the probability of a specific individual located in C to be chosen as the ancestor in a birth event is slightly elevated. This spatial influence is of course not present in the BD.

We simulated 100 runs of the ΛV forward in time with the choices for θ and ξ as in the previous section such that θ 2 ξ = 1, and again υ = 1. We still considered a rectangle of size 10 × 10 and the population density was set to ρ = 4. Under these assumptions, at time 0 the number of individuals on the rectangle is Poisson-distributed with mean 400 (i.e., ρwh). After T = 1 and T = 4 units of simulated time, we recorded the distribution of family sizes and formed the average over all runs. The same was done for a BD with λ = µ = 2π.

The frequencies of family sizes under the BD generally agree well with Eq. ( 9). Hence we represent the BD by this function in Figure 2. For the ΛV process, the absolute values of θ and ξ visibly affect the shape of the distribution. If θ is large and events comparably rare (e.g., in the setting θ = 1, ξ = 1), we observe an overabundance of families of size one, and a much flatter distribution otherwise, with extended frequencies of higher family sizes. This observation is most likely explained by the variance in offspring number when a REX event takes place, which is given by 2πθ 2 ρυ and thus quadratic with respect to θ. Smaller values of θ typically provide a good fit between ΛV and BD. However, after four units of time, we observe a deficit in large family sizes in ΛV versus BD. That discrepancy probably reflects the impact of spatial constraints in the ΛV. Indeed, families with most members located close to a boundary give birth to a smaller number of individuals compared to those located far away from these boundaries. This difference of behaviour is probably responsible, at least in part, for the observed divergence between the two models although additional investigations are clearly needed in order to have a deeper understanding of the forces at play.

Properties of BD and ΛV * as tree-generating processes

The statistical properties the BD model as a tree-generating process are well-known, e.g., with respect to branch lengths and tree topology (see e.g., [START_REF] Gernhard | Stochastic models for speciation events in phylogenetic trees[END_REF]). In particular, consider the following setting: Assume that the BD process is initialised at t or units of time in the past with a single lineage, and there are n > 0 lineages alive at the present (time t 0 = 0). Consider the joint density p BD (τ, t 1 , . . . , t n-1 | n, λ, µ, t 0 , t or ) [START_REF] Guindon | Demographic inference under the coalescent in a spatial continuum[END_REF] of the topology τ and the bifurcation times T 1 , . . . , T n-1 in the past of the family genealogy (where T 1 is the most recent bifurcation and T i < T i+1 ), given the family size n, the time frame [0, t or ] and the parameters of the process. Then, it holds that

p BD (τ, t 1 , . . . , t n-1 | n, λ, µ, t 0 , t or ) ∝ n-1 i=1 p 1 (t i ) (11)
where p 1 (t) is the probability that a BD process starting at time 0 with one lineage has again one single lineage after t units of time (see e.g., [START_REF] Yang | Bayesian phylogenetic inference using dna sequences: a markov chain monte carlo method[END_REF]). For a critical BD process (λ = µ), we have

p 1 (t) = 1 (1 + λt) 2 (12)
We compare p BD (τ, t 1 , . . . , t n-1 | λ, µ, t 0 , t or , n) to p ΛV (τ, t 1 , . . . , t n-1 | ξ, θ, t 0 , t or , n) through simulations in the case where n = 2. We generated trees under the ΛV process forward in time using the following procedure: the value of t or is chosen arbitrarily and the corresponding initial location is chosen uniformly at random in A. We run the process, updating the genealogy of descendants of the founder after each REX event, until time 0 is reached. Simulations are discarded whenever the number of lineages n is different from two. We retain a sample of genealogies with valid realisations of T 1 . We then compared the empirical distribution of this random variable to that derived analytically for the BD. We repeated these simulations for different values of θ, with ξ chosen such that ξ = 1/θ 2 , and therefore λ = 2π. We opted for t or = 0.5, since this suffices to outline the shape of T 1 for the range of values of θ selected here.

Figure 3 shows that for large radii (θ = 1, in particular), the distribution of coalescence times of two lineages noticeably diverge in shape and mode from that derived from the BD process. We hypothesise that the number of REX events involved in these particular simulation settings is relatively small so that lineages have to "wait" relatively long periods of time before being affected by an event, preventing early coalescent events. For smaller values of θ (and therefore larger values of ξ), distributions of T 1 derived from the ΛV are more similar to that given by the BD, as expected.

We now focus on n = 3 and compare the bifurcation times T 1 and T 2 obtained under the ΛV and the BD processes. For the joint density of the split times in a critical BD conditioned on n = 3 and t or , it holds that

p BD (τ, t 1 , t 2 | n = 3, λ, µ, t 0 , t or ) ∝ 1 (1 + λt 1 ) 2 • 1 (1 + λt 2 ) 2 (13)
Since t 1 < t 2 , we can calculate the marginal density p BD (t 1 | n, λ, t 0 , t b ) for T 1 under the above conditions and using [START_REF] Hastings | Monte carlo sampling methods using markov chains and their applications[END_REF]: The distributions for the ΛV were obtained from simulations with 100 repeats for each value of θ. The densities corresponding to the BD (in black) agree with Equations ( 14) and ( 16).

p BD (t 1 | n, λ, t 0 , t or ) ∝ p 1 (t 1 ) tor t 1 p 1 (t 2 )dt 2 (14) = 1 (1 + λt 1 ) 2 • 1 λ + λ 2 t 1 - 1 λ + λ 2 t or (15)
Similarly, we obtain the marginal density p BD (t 2 |n, λ, t or ) of t 2 :

p BD (t 2 | n, λ, t 0 , t or ) ∝ p 1 (t 2 ) t 2 0 p 1 (t 1 )dt 1 (16) = 1 (1 + λt 2 ) 2 • 1 λ - 1 λ + λ 2 t 2 (17)
As for the ΛV, we repeated the simulations described above, this time discarding all instances where the number of lineages n was not equal to three at time 0. Also, we discarded cases where two of the three final lineages were generated in the same birth event as in such a case τ is not a binary tree. However, with decreasing θ, this typre of event becomes less and less likely.

The bifurcation times t 1 and t 2 can then given by the genealogies obtained from the successful runs. We used the same parameter combinations as in the case of n = 2, except for θ = 0.01 and ξ = 10, 000, as according to our observations, this case becomes numerically infeasible to simulate in a reasonable amount of computing time. Here, the starting point of the simulations was taken as t or = 2 units of time in the past.

Results in Figure 4 indicate a good agreement between distributions of T 1 and that of T 2 for the two models for values of θ smaller than 1. Although obtaining a sufficiently large number of valid draws from the target distributions was computationally challenging (hence the rough aspect of some of the curves derived from ΛV simulations), the modes of the reconstructed densities get closer to that of the BD process when the radius decreases.

Comparison of ΛV and BD 2 processes

Results in the previous section indicate that the ΛV and the BD treegenerating processes are, at least in the simulation settings examined in the present study, equivalent in the limit of a small radius and a large rate of REX events. The present section aims at assessing whether the similarity between the two models still stands when incorporating spatial information.

When considering a single lineage and ignoring border effects, the movements of the corresponding particle evolving under ΛV follows a (shifted) Brownian process with diffusion parameter 4πθ 4 ξ (see Appendix, section 4.1). The behaviour of a pair of lineages is not as straightforward as that of two independent Brownian trajectories. In particular, during the period of time following the birth of the two lineages (i.e., moments after the splitting of their ancestor), the two particles remain in the vicinity of one another. Any given event affecting one of the two particles is thus likely to impact the other as well. The movements of the two particles are therefore not independent and the correlation depends on the time to their common ancestor. Yet, in the limit of a small radius, one may expect the dependency between particles to vanish quickly after their birth and particles may thus be considered as independent when monitored over relatively long periods of time. However, the impact of borders in the habitat can no longer be ignored under the ΛV while these do not play a role in the BD 2 . The next sections explore these issues using forward and backward in time simulations of two lineages under both processes.

Comparison of likelihoods

We first focus on the comparison of both models by considering their respective predictions of the spatial coordinates at the tips of a two-lineage tree with fixed ancestral node age and location. The density of interest is noted here as q f (L 2 , L 3 | t 1 , l 1 , θ, ξ, υ, h, w, n = 2), corresponding to the joint density of coordinates L 2 and L 3 at the tips of lineages 2 and 3, given the time t 1 at which these two lineages coalesce, l 1 the location of the ancestor just before the edge splitting event and the parameters of the ΛV model (with λ = µ = 2π and σ 2 = 4πθ 2 c for BD 2 ). The subscript f in the density stands for "forward in time".

The habitat is modelled as a 10 × 10 square (i.e., h = w = 10) with an ancestral location set to l 1 = (5, 5). The radius θ is fixed to 0.025 throughout these simulations and the rate of events ξ is equal to 1/θ 2 so that c = 1, as For each model, we generated 1,000 draws from the corresponding distribution with density q f (L 2 , L 3 | t 1 , l 1 , θ, ξ, υ, h, w, n = 2), for values of t 1 = 100 (left), 1,000 (centre) and 5,000 (right), with θ = 0.025 and ξ = 1/θ 2 . The density plots display the joint distributions of L 2 and L 3 along the x-axis (denoted l 2 (x) and l 3 (x) respectively). Both axes have lower and upper limits -5.0 and +15, to be compared with limits of the habitat (i.e., lower and upper limits of 0 and 10 along both axes) per usual. We then obtained the joint distributions of L 2 and L 3 for values of t 1 equal to 100, 1,000 and 5,000.

Figure 5 shows that for relatively small values of the coalescence time, both models predict virtually identical distributions of locations at the tips. In other words, tip locations under the ΛV are well approximated by a multivariate normal when the radius of events is small compared to the size of the habitat and the rate of REX events is large. For larger values of t 1 , border effects impact ΛV substantially and the BD 2 model puts a large probability mass on tip locations falling outside the habitat (see Figure 5 right). In these conditions, the distribution of L 2 and L 3 under ΛV becomes almost uniform and is thus clearly distinct from a bivariate normal (even in the case where realisations of L 2 and L 3 under BD 2 are generated using a bivariate normal truncated to [0, h] × [0, w] so as to better accommodate for the limits of the habitat (results not shown)).

Comparison of posterior densities

Results obtained in section 2.4.1 indicate that the likelihood of both models are only equivalent in cases where the time to coalescence is not too distant in the past so that the impact of the limits of the habitat can be safely ignored. We now focus on the distribution of the coalescence time and the corresponding ancestral location conditioned on the sampled locations of the two focal lineages. Let q b (L 1 , T 1 | l 2 , l 3 , θ, ξ, υ) denote that distribution, with the subscript b for the "backward in time" process. The forward (q f (•), see previous section) and backward (q b (•)) distributions bare obvious connections (see below). Yet, just because recent coalesence times most likely generate pairs of tips that are located in a small area (see Figure 5 left) does not necessarily imply that the most probable times of coalescence of lineages sampled in such region are young.

We generated samples from the target distribution through direct simulation under the ΛV model (see section 2.5.2). As for the BD 2 model, we obtained correlated samples by applying a Metropolis-Hastings algorithm [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte carlo sampling methods using markov chains and their applications[END_REF] with standard proposal operators for updating the time to coalescence and the corresponding spatial coordinates. Figure 6 shows the distribution function of T 1 and L 1 (focusing on the x-axis) obtained under the two models for tip coordinates set to l 2 = (5.00, 5.43) and l 3 = (4.75, 5.00), and habitat size defined using w = h = 10. We considered a similar range of values for the radius as the one used previously, i.e., θ = 1, 0.25, 0.05 and 0.025. The distribution of T 1 shows a behaviour similar to that observed when ignoring spatial information with cumulative distributions of the two models becoming more similar as the radius decreases. Results obtained for the spatial component of the models are noticeably different. The range of values for L 1 is much narrower under BD 2 compared to ΛV, with a distribution converging to coordinates generally tightly grouped midway between the two sampled tip locations (i.e., (5.00 + 4.75)/2 along the x-axis), while the ΛV shows a much broader distribution of estimated ancestral locations, even though an inflexion of the distribution function is observed as well around the midpoint between the sampled lineages.

At first glance, the comparison of results obtained by running BD 2 forward and backward may appear puzzling: forward in time simulations show considerable variance in tip coordinates (see Figure 5, right) while backward in time simulations, starting from the most likely tip locations and considering time to coalescent of the same order of magnitude, yields very precise coordinates at the coalescent node (see Figure 6, right). For a fixed coalescent time t 1 , the posterior distribution of the spatial coordinates at the coalescent node is derived as follows:

q b (l 1 | l 2 , l 3 , t 1 , σ) ∝ q f (l 2 , l 3 | l 1 , t 1 , σ) ∝ ϕ(l 2 ; l 1 , σ 2 t 1 )ϕ(l 3 ; l 1 , σ 2 t 1 ) = ϕ(l 1 ; l 2 + l 3 2 , σ 2 t 1 2 )
where ϕ(•, µ, σ 2 ) is the normal density with mean µ and variance σ 2 . Although the following statement lacks a sound mathematical backing, one may argue that the diffusion parameter of the backward in time process is thus half that of the forward in time process. This observation explains, at least partially, the difference of behaviour of the forward and backward versions of the BD 2 . Explaining the differences between the BD 2 and the ΛV models is less straightforward. Conditioned on the time to coalescence of the two focal lineages, the distribution of REX events is no longer uniform in space, prohibiting simple mathematical results about the spatial coordinates on the ancestor. The simulation results presented in this study simply suggest that, when focusing on the spatial component of the models, the ΛV and BD 2 behave differently in the limit of small radius and frequent REX events even though both models are equivalent when focusing on a single lineage.

Efficient simulations under the SLFV model

Comparison between the ΛV and other tree-and spatial coordinates-generating processes depends on our ability to efficiently simulate data under these stochastic models. In particular, for the ΛV model, the present study required simulation under the backward and forward in time versions of this process. As seen above, the forward process generates realisations that may be used for direct comparison with the likelihood of the model. Backward simulations are used instead for comparison with the posterior densities of ancestral node ages and their spatial coordinates. Since we focus on the limit of small radius in the present study, the vast majority of REX events do not hit any lineage, making the simulations computationally inefficient. For instance, the backward generation of a two-lineage data set with θ = 0.02 takes about 45 minutes for lineages that are 0.5 space unity away from each other on a 10 × 10 square. Also, naive forward in time simulations require to monitor the whole population of lineages and keep track of their positions at each REX event, which is costly in terms of memory usage. We provide below two algorithms for forward and backward simulation of two lineages evolving under the ΛV process that alleviate these difficulties.

Forward simulations

Our objective here is to obtain independent random draws from the distribution with density q f (L 2 , L 3 | t 1 , l 1 , θ, ξ, υ, n = 2). In words, we want to generate locations L 2 and L 3 for two focal lineages (2 and 3, sampled at time 0) given that their most recent common ancestor split at time t 1 and had location l 1 just before the split. In the following, we first give an algorithm that simulates the trajectory of two lineages forward in time which does not require monitoring the whole population. We then describe a modified, more efficient version of this method that ignores events that leave the two lineages unchanged.

We first generate the position z of the REX event corresponding to the split of the lineage ancestral to 2 and 3 by sampling from a truncated normal distribution with mean l 1 , variance θ 2 and truncation set so that z falls within the h × w rectangle defining the habitat. Next, we choose the initial position of each of the two focal lineages, noted l 2 and l 3 , by sampling from a truncated normal with mean z and variance θ 2 . (1) The time to the next REX event is then obtained by sampling from an exponential distribution with rate ξwh. (2) The position of that event is selected uniformly at random in the h×w rectangle. (3) The probability that the sampled lineage i is hit by this event is

u i (z) = υ exp -||l i -z|| 2 θ 2
(noted as u i in the following), where i = 2 or 3. Also, the probability that both lineages are hit is u * = u 1 u 2 . We need to exclude the situation where both lineages are hit by the REX event, since only one parent is selected to give birth to new lineages per REX event. If both were hit by one event, at least one of the two lineages would die without offspring and would therefore not survive to the present time (t = 0). ( 4) The probability that one and only one of the two lineages is hit is thus u 1 + u 2 -2u * . If this event takes place, it affects i with probability u i -u * u 1 +u 2 -2u * and the new position of lineage i is sampled from a truncated normal with mean c and variance θ. Steps (1)-(4) of the above procedure are repeated until the time elapsed, i.e. the sum of exponentially distributed times generated in (1), exceeds t 1 .

The present study focuses on the case where the radius of events is small compared to the size of the habitat. As already mentioned, in this situation, most events do not impact any of the sampled lineages, conveying limited information for our purpose (the rate of these events enters the model as a time scaling factor). We thus elected to adapt our simulation procedure so as to focus solely on the rate of events where one sampled lineage and only one dies. This rate is simply the product of the rate of all events (ξhw) by the probability that one of the two lineages dies, i.e., 1 hw z∈A (u 1 + u 2 -2u * )dz, where A is the h×w rectangle and therefore varies with the lineages' positions (see Appendix for the solution to that integral).

When focusing only on events that impact the sampled lineages, the spatial position of the event centres is no longer uniform. Deriving the joint distribution of the REX centre position along with that of the two lineages right after the event is thus essential in designing an approach that generates random draws from the correct distribution. Although the ordering in which lineages are considered when examining the impact of an event is not relevant, we hereby consider our two focal lineages in a serial fashion, i.e., one lineage is considered as the first while the other is the second. Let H 1 be a discrete random variable with state space {2, 3} corresponding to the event space {"lineage 2 is the first lineage and dies","lineage 3 is the first lineage and dies"}. Also, let (H 2 | z) be the random variable with state space {1, 2} corresponding to the event space {"the second lineage dies", "the second lineage does not die"}. The probability density of interest is thus noted as:

p(H 1 = 2, H 2 = 2, z | θ) + p(H 1 = 3, H 2 = 2, z | θ) = Pr(H 1 = 2)p(z | H 1 = 2, θ) Pr(H 2 = 2 | z, H 1 = 2)+ Pr(H 1 = 3)p(z | H 1 = 3, θ) Pr(H 2 = 2 | z, H 1 = 3) = 1 2 p(z | H 1 = 2, θ) Pr(H 2 = 2 | z)+ 1 2 p(z | H 1 = 3, θ) Pr(H 2 = 2 | z)
Examination of the last expression suggests that the following procedure could be used in order to get a valid random draw for the event centre according to the model of interest: (1) pick one of the two lineages uniformly at random as the first lineage. Let i denote the event corresponding to the death of that lineage; (2) sample the value of z from the distribution with density p(z | H 1 = i, θ); (3) let u be a random draw from

U [0, 1], if u ≤ Pr(H 2 = 2 | z, H 1 = i) (i.
e., the second lineage dies), return to (1), otherwise return z.

Backward simulations

The goal of the backward simulations is to generate independent random draws from the distribution with density q(T 1 , L 1 | l 2 , l 3 , θ, ξ, υ), i.e., given l 2 and l 3 , the locations of the two sampled lineages at present, generate T 1 and L 1 , the time and location of their most recent common ancestor. As noted above, if done naively, simulation of the ΛV when tracking a small number of sampled lineage is computationally costly since the vast majority of REX events do not impact any of the sampled lineage. A more efficient approach would then be to focus exclusively on the REX events that either hit one lineage only, or hit both of them as is the case when coalescence take place, and set the rate of these events in an appropriate manner. Below is a description of one such approach.

The rate of events that hit one or the two lineages is given by the product of the rate of all types of events (ξwh) by the probability that one or the two lineages are hit, i.e., using the notation from the previous section:

1 wh z∈A (u 1 + u 2 -u 1 u 2 )dz (see section 4.2 in the Appendix). Hence, here again, this rate is not constant in time as it changes with the position of lineages. The core of the proposed procedure relies on the distribution of the location of a REX event conditioned on that event hitting both lineages or only one of them. Using a similar approach as for the forward case, let H 1 be a discrete random variable with state space {2, 3} corresponding to the event space {"lineage 2 is the first lineage and is hit", "lineage 3 is the first lineage and is hit"}. Also, let (H 2 | z) be the random variable with state space {1} corresponding to the event space {"the second lineage is hit or not"}. The joint probability density of one or the two lineages being hit by the event and the location of the REX event is thus expressed as follows:

p(H 1 = 2, H 2 = 1,z | θ) + p(H 1 = 3, H 2 = 1, z | θ) = Pr(H 1 = 2) × p(z | H 1 = 2, θ) × Pr(H 2 = 1|z, H 1 = 2)+ Pr(H 1 = 3) × p(z | H 1 = 3, θ) × Pr(H 2 = 1|z, H 1 = 3) = Pr(H 1 = 2) × p(z | H 1 = 2, θ) + Pr(H 1 = 3) × p(z | H 1 = 3, θ)
The last expression above suggests that simulating a valid value for the centre position can be done by first picking one of the lineages to be hit by the event with probability Pr(H 1 = •) = 1/2 and then sampling the event centre from a truncated normal centred on that lineage (with variance θ 2 ), i.e., with the corresponding density p(z | H 1 = •, θ). The simulation continues if the second lineage is not hit by the same event. It stops if the second lineage is hit by the event. In the second case, one then samples L 1 from a truncated normal centred on z and the simulation is complete.

Discussion

The present study illustrates several parallels between the ΛV and BD models. Starting from the observation that in the ΛV lineages experience birth and death events over the course of time in a manner similar to the BD, we derived analytical results concerning the rates of these events in the ΛV when approaching the limit θ → 0, ξ → ∞ and θ 2 ξ → c for constant c (we use ΛV* to denote this particular version of the ΛV). We verified through simulations the theoretical predictions and investigated several related questions regarding the genealogical process in the ΛV. The ΛV was simulated backward and forward in time in accordance with its standard formulation [START_REF] Barton | A new model for evolution in a spatial continuum[END_REF][START_REF] Barton | Modelling evolution in a spatial continuum[END_REF]. We introduced two algorithms that permit efficient simulation by skipping REX events that do not impact the sampled lineages. We also implemented forward and backward numerical techniques, through direct simulation or the sampling of correlated samples through MCMC, under the BD tree-generating process and BD with Brownian evolution of spatial coordinates along the tree edges (the so-called BD 2 model).

We first focused on the tree generating process induced by the ΛV* model. Our simulations indicate that the per lineage birth and death rates do indeed converge to that derived analytically, thereby establishing a first connection with the BD model. We next focused on the distribution of the number of descendants of individual lineages after fixed amounts of time. Here again, we observe a good agreement between the two models, especially for short waiting times. The distributions become distinct for longer time periods, at which point the size of surviving families is large so that the effect of the limited size of the habitat cannot be ignored under the ΛV model, while it plays no role under the BD model. Finally, forward simulations suggest that the times to first and second coalescent events in samples of size three in the ΛV converge in distribution to those observed in the BD. Altogether, our results indicate that the tree-generating processes induced by the ΛV* and BD processes are equivalent as long as the sample size is small enough so that the limits of the habitat can safely be ignored.

When spatial coordinates of lineages are taken into account, the finite rectangle we simulate on with the ΛV process induces boundary effects, causing a differentiation between the densities of ancestral lineage locations for the ΛV and the BD 2 models. This discrepancy does not vanish with larger rates of REX events and smaller radius. Here, the impact of a decreasing radius does not seem to be offset by the increasing rate of events, pushing coalescent times deeper in the past, thus making the probability for any lineage to hit the habitat boundaries before coalescing non negligible. Backward in time simulations of the dynamics of a pair of lineages show that the spatial distribution of the most recent common ancestor is substantially less variable under BD 2 compared to ΛV*. This observation entails serious consequences in practice as it implies that the choice of model will impact on the precision with which ancestral coordinates are to be estimated, with BD 2 potentially giving overly precise estimates when the model that actually generated the data of interest is closer to ΛV.

Finally, we present two algorithms for simulating the temporal and spatial dynamics of a pair of lineages forward and backward in time. These new methods are computationally efficient as they focus solely on REX events that impact the lineages under scrutiny while the naive approach simulates vast numbers of events affecting individuals in the population that are not incorporated in the sample. Importantly, the new backward in time algorithm may serve as a basis for the simulation-based inference of model parameters under the ΛV (using, for instance, approximate Bayesian computation). While the ΛV model is amenable to parameter inference [START_REF] Guindon | Demographic inference under the coalescent in a spatial continuum[END_REF], the task is computationally challenging. Efficient approximation for the time to coalescence of pairs of lineages were derived recently [START_REF] Wirtz | Rate of coalescence of lineage pairs in the spatial λ-Fleming-Viot process[END_REF]. Yet, fast and accurate parameter estimation methods are still lacking and the proposed simulation algorithm presented in this study may contribute to filling this void.

= πυθ 2 2|A| erf

√ 2z x 2θ -erf √ 2(z x -w) 2θ × erf √ 2z y 2θ -erf √ 2(z y -h) 2θ . ( 19
)
In cases where the argument of each error function above is large enough (i.e., greater than ≃ 2), its value is close to one. These conditions are met when θ ≪ min(z x , z y ) and z is far enough from the edges of the habitat (i.e., w -z x ≫ 0 and z y ≫ 0, and likewise for z y ) . In this situation, the expression above simplifies, yielding

u(l + , z) |A| d 2 l + ≃ 2πυθ 2 /|A|, (20)
which is also the marginal probability of the lineage being hit (i.e., without conditioning on the position of the REX event). We will consider that this approximation holds in what follows. The rate at which a given lineage is hit is thus 2ξπθ 2 υ.

Also, the probability density of l -(the position of the lineage just after the REX event, still going backward in time) given l + (with l -̸ = l + ) is 1 4π 2 θ 4 v(l -, z i )v(l + , z i )d 2 z i . This integral yields 1 4πθ 2 exp -1 4θ 2 ||l --l + || 2 , i.e., a bivariate normal density with mean l + and covariance matrix 2θ 2 I. The variance of offspring location in a one-dimensional space given the parental location is thus E(d 2

x ) = 2θ 2 ; where d 2 x is the squared Euclidean distance in a one dimensional habitat. θ 2 is thus half the expected square Euclidean distance between parent and offspring in one dimension. In two dimensions, we have E 1 2 (d 2 x + d 2 y ) = 1 2 E(d 2 x ) + E(d 2 y ) = 2θ 2 i.e., θ 2 is a quarter of the expected square Euclidean distance between parent and offspring. In a n-dimensional habitat, θ 2 is 1/2n times this expected distance.

Altogether, in a two-dimensional habitat, the variance of spatial coordinates of a lineage along a given axis thus increases with time proportionally to σ 2 := 4θ 4 ξπυ. In the limit where λ → ∞ and θ → 0, we hypothesise that the backward-in-time motion of a single lineage is a Brownian process with diffusion parameter σ 2 .

Probability of coalescence of two lineages

Let l 2 = (l 2,x , l 2,y ) and l 3 = (l 3,x , l 3,y ) be the current positions of the two lineages under scrutiny. The probability that lineage i is hit given the centre position c = (c x , c y ) is, by definition of the ΛV model, u i (c) = υ exp -||l i -c|| 2 2θ 2 . The probability that both lineages are hit (i.e., coalesce) is obtained following an approach similar to that used for a single lineage (see above): for small values of θ.

Figure 1 :

 1 Figure 1: Densities of the times until a given individual is subject to a birth event (left) and to a death event in the ΛV. Tested values of θ are given to the right. ξ is such that θ 2 ξ = 1. The respective densities in the BD are shown in black.

Figure 2 :

 2 Figure 2: Distribution of family sizes under BD and ΛV for T = 1 (left column) and T = 4 (right column) for various radii. The distribution were obtained analytically for BD (see main text) and simulated forward in time for ΛV.

Figure 3 :Figure 4 :

 34 Figure3: Distributions of T 1 for two-tip trees under the ΛV and the BD processes. The distributions for ΛV were obtained from simulations with 100 repeats for each value of θ while that for the BD (in black) is analytical (see main text).

Figure 5 :

 5 Figure5: Distributions of tip locations under the ΛV and BD 2 models. For each model, we generated 1,000 draws from the corresponding distribution with densityq f (L 2 , L 3 | t 1 , l 1 , θ, ξ, υ, h, w, n = 2), for values of t 1 = 100 (left), 1,000 (centre) and 5,000 (right), with θ = 0.025 and ξ = 1/θ 2 . The density plots display the joint distributions of L 2 and L 3 along the x-axis (denoted l 2 (x) and l 3 (x) respectively). Both axes have lower and upper limits -5.0 and +15, to be compared with limits of the habitat (i.e., lower and upper limits of 0 and 10 along both axes)

Figure 6 :

 6 Figure6: Posterior distribution functions of coalescent times (top row) and spatial coordinates (bottom row) under the ΛV and BD 2 models. Samples from the distribution with density q b (L 1 , T 1 | l 2 , l 3 , θ, ξ, υ) were generated under both models for values of θ =1, 0.25, 0.05 and 0.025 (with ξ = 1/θ 2 and υ = 1). Tip coordinates for the two sampled lineages were set to (5.00,5.43) and (4.75, 5.00).

  on the REX event having location z = (z x , z y ) isu(l + , z) |A| d 2 l + = + , z) |A| dl x dl y

Pr

  x -z x ) 2 + (l 2,y -z y ) 2 2θ 2 × exp -(l 3,x -z x ) 2 + (l 3,y -z y ) 2 2θ 2 dz x dz y = πθ 2 υ 2 4|A| exp -(l 2,x -l 3,x ) 2 + (l 2,y -l 3,y ) 2 4θ 2 × erf l 2,x + l 3,x 2θ -erf l 2,x + l 3,x -2w 2θ erf l 2,y + l 3,y 2θ -erf l 2,y + l 3,y -2h 2θand the probability of coalescence gets close to πθ 2 υ 2 |A| exp -||l 2 -l 3 || 2 4θ 2

Acknowledgements

This work was financially supported by the Agence Nationale pour la Recherche [https://anr.fr/] through the grant GENOSPACE, and the Walter-Benjamin Program (WI 5589/1-1) of the DFG [https://dfg.de/].

Appendix

Dispersal of a single lineage under ΛV

When considering the backward in time ΛV process, the rate at which a lineage is hit by a REX is the product of the rate at which these events occur (ξwh = ξ|A|) by the probability that a lineage is hit. Let l + be the (twodimensional vector) location of the focal lineage just before the REX event that occurred at time t. The probability that this lineage is hit conditional


Mona Alleshagh

Samahe Sadjadi
email: s.sadjadi@ippi.ac.ir

Hassan Arabi

Naeimeh Bahri-Laleh

Eric Monflier

Pd on ligand decorated chitosan as an efficient catalyst for hydrofinishing of polyalphaolefin: combination of experimental and computational study

Keywords: Chitosan, Hydrofinishing, Polyalphaolefins, DFT

As hydrofinishing of polyalphaplefins is a crucial step for improving their features and expanding their industrial applications, in this research a novel chitosan-based catalyst was designed for this purpose. To synthesize the catalyst, chitosan was first modified with a ligand. In this regard, computational study has been employed for devising an efficient ligand from a series of possible candidates. Then, chitosan was functionalized with the selected ligand, palladated and used as a heterogeneous catalyst for hydrofinishing of polyalphaplefin. The reaction variables, including reaction temperature, pressure and catalyst loading were optimized and it was found that 5 wt.% catalyst could promote the hydrofinishing reaction under hydrogen pressure of 6 bar at 130 ºC to give the hydrogenated oil in 98% yield. Moreover, the nature of catalysis was heterogeneous and the catalyst could be recovered and reused for five runs with slight loss of its catalytic activity and Pd leaching.

Introduction

Chitosan, CS, is a natural polymer [START_REF] Hajipour | Magnetic chitosan-functionalized cobalt-NHC: Synthesis, characterization and catalytic activity toward Suzuki and Sonogashira cross-coupling reactions of aryl chlorides[END_REF][START_REF] Jiao | Recyclable copper catalyst on chitosan for facile preparation of alkyl/aryl mixed phosphates via deaminated esterification between diphenylphosphoryl azides and aliphatic alcohols[END_REF][START_REF] Anbu | Knoevenagel-Doebner condensation promoted by chitosan as a reusable solid base catalyst[END_REF][START_REF] Maleki | Green multicomponent synthesis of four different classes of six-membered N-containing and O-containing heterocycles catalyzed by an efficient chitosan-based magnetic bionanocomposite[END_REF] that has vast range of applications [START_REF] Sahu | Chitosan: An Efficient, Reusable, and Biodegradable Catalyst for Green Synthesis of Heterocycles[END_REF][START_REF] Sahu | Correction for "Chitosan: An Efficient, Reusable, and Biodegradable Catalyst for Green Synthesis of Heterocycles[END_REF][START_REF] Vincent | Chitosan-Supported Palladium Catalyst. 1. Synthesis Procedure[END_REF] in various research domains [START_REF] Torkaman | Modification of chitosan using amino acids for wound healing purposes: A review[END_REF][START_REF] Araujo | Chitosan-based systems aimed at local application for vaginal infections[END_REF][START_REF] Cavallaro | Chitosan-based smart hybrid materials: a physico-chemical perspective[END_REF][START_REF] Baran | Pd/CoFe2O4/chitosan: A highly effective and easily recoverable hybrid nanocatalyst for synthesis of benzonitriles and reduction of 2-nitroaniline[END_REF][START_REF] Hussein | Effect of TiO2 nanoparticles on the structural configurations and thermal, mechanical, and optical properties of chitosan/TiO2 nanoparticle composites[END_REF]. This carbohydrate that can be achieved through deacetylation of chitin [START_REF]Heterogeneous catalysis on chitosan-based materials: a review[END_REF][START_REF] Bodhak | An efficient and recyclable chitosan supported copper (II) heterogeneous catalyst for C-N cross coupling between aryl halides and aliphatic diamines[END_REF][START_REF] Shen | A highly active and easily recoverable chitosan@ copper catalyst for the C-S coupling and its application in the synthesis of zolimidine[END_REF], has been extensively applied for the catalytic purposes [START_REF] Anbu | Chitosan as a biodegradable heterogeneous catalyst for Knoevenagel condensation between benzaldehydes and cyanoacetamide[END_REF][START_REF] Koohestani | Polyionic liquid decorated chitosan beads as versatile metal-free catalysts for catalyzing chemical reactions in aqueous media[END_REF][START_REF] Macquarrie | Applications of Functionalized Chitosan in Catalysis[END_REF]. The main advantage of CS is the presence of ample functionalities, both -OH and amine groups, in its backbone [START_REF] Sadjadi | Ionic liquid decorated chitosan hybridized with clay: A novel support for immobilizing Pd nanoparticles[END_REF][START_REF] Wu | EDTA modified βcyclodextrin/chitosan for rapid removal of Pb (II) and acid red from aqueous solution[END_REF] that allows functionalization of CS and incorporation of further functionalities [START_REF] Zeng | N-doped mesoporous carbons supported palladium catalysts prepared from chitosan/silica/palladium gel beads[END_REF][START_REF] Zeng | Novel macroporous palladium cation crosslinked chitosan membranes for heterogeneous catalysis application[END_REF][START_REF] Dhanavel | Preparation and characterization of cross-linked chitosan/palladium nanocomposites for catalytic and antibacterial activity[END_REF][START_REF] Xie | A chitosan modified Pt/SiO2 catalyst for the synthesis of 3-poly (ethylene glycol) propyl ether-heptamethyltrisiloxane applied as agricultural synergistic agent[END_REF][START_REF] Akpomie | Mechanism on the sorption of heavy metals from binary-solution by a low cost montmorillonite and its desorption potential[END_REF]. On the other hand, Similar to many supports [START_REF] Maleki | Pumice-modified cellulose fiber: An environmentally benign solid state hybrid catalytic system for the synthesis of 2,4,5triarylimidazole derivatives[END_REF][START_REF] Maleki | Green oxidation protocol: Selective conversions of alcohols and alkenes to aldehydes, ketones and epoxides by using a new multiwall carbon nanotube-based hybrid nanocatalyst via ultrasound irradiation[END_REF], multiple heteroatoms in CS structure provide electrostatic interactions with other species, such as metallic nanoparticles (NPs). This characteristic can be exploited for efficient stabilization of catalytic species and development of bio-based heterogeneous catalysts [START_REF] Zeng | Microstructure-stability relations studies of porous chitosan microspheres supported palladium catalysts[END_REF][START_REF] Zeng | Chitosan microspheres supported palladium heterogeneous catalysts modified with pearl shell powders[END_REF][START_REF] Zeng | Novel chitosan-based/montmorillonite/palladium hybrid microspheres as heterogeneous catalyst for Sonogashira reactions[END_REF][START_REF] Zhao | Chitosan supported Pd0 nanoparticles encaged in Al or Al-Fe pillared montmorillonite and their catalytic activities in Sonogashira coupling reactions[END_REF][START_REF] Leonhardt | Chitosan as a support for heterogeneous Pd catalysts in liquid phase catalysis[END_REF][START_REF] Shuo | Applications of Chitosan and Its Derivatives Supported Transitional Metal Catalysts in Organic Reactions[END_REF].

Low-viscosity lubricants can impressively reduce friction loss in engines, thereby making them energy efficient, which will impart to global energy consumption saving. One of the major lubricant families of synthetic oils type, used to reach the so-called goal is polyalphaolefins (PAO type) lubricants. PAOs belong to the group IV of American Petroleum Institute (API) category.

They are hydrogenated α-olefin (primarily C8-C12) oligomers [START_REF] Gharajedaghi | Experimental and DFT study on titanium-based half-sandwich metallocene catalysts and their application for production of 1-hexene from ethylene[END_REF][START_REF] Hanifpour | 1-Decene oligomerization by new complexes bearing diamine-diphenolates ligands: Effect of ligand structure[END_REF][START_REF] Hanifpour | Group IV diamine bis(phenolate) catalysts for 1-decene oligomerization[END_REF], used as base oil in highperformance engine lubricants [START_REF] Hanifpour | Coordinative chain transfer polymerization of 1-decene in the presence of a Ti-based diamine bis(phenolate) catalyst: a sustainable approach to produce low viscosity PAOs[END_REF]. In fact, PAOs have all hydrocarbon structure, without double bonds, ring structures, S and N containing components or waxy hydrocarbons [START_REF] Fallah | Interaction of common cocatalysts in Ziegler-Natta-catalyzed olefin polymerization[END_REF]. The non-polar nature of PAOs, together with their high viscosity index (VI, >125), good oxidation stability, low pour point and good compatibility with mineral oils make them prominent candidates in engine oil topic, at both low and high temperature applications [START_REF] Zhang | Preparation, characterization and tribological properties of polyalphaolefin with magnetic reduced graphene oxide/Fe3O4[END_REF].

The current industrial process to get final PAO product, involves a 3-stages process consisting of oligomerization, hydrogenation and distillation. In fact, crude PAOs bear C=C double bonds, which are susceptible to the oxidation at high temperature of usage [START_REF] Hanifpour | Kinetic and microstructural studies of Cp2ZrCl2 and Cp2HfCl2-catalyzed oligomerization of higher α-olefins in mPAO oil base stocks production[END_REF]. To overcome this drawback, hydrofinishing of crude PAO is accomplished, usually at high pressures and temperatures (> 30 bar and >200 °C). This extremely increases the risk and cost of PAO production. Therefore, industry is always looking for developed and efficient catalytic systems operate at mild reaction conditions in the hydrogenation process. One of the mostly utilized catalyst is Pd catalyst [START_REF] Bahri-Laleh | Pd immobilized on dendrimer decorated halloysite clay: Computational and experimental study on the effect of dendrimer generation, Pd valance and incorporation of terminal functionality on the catalytic activity[END_REF][START_REF] Maleki | Palladium-decorated ophenylenediamine-functionalized Fe3O4/SiO2 magnetic nanoparticles: A promising solid-state catalytic system used for Suzuki-Miyaura coupling reactions[END_REF][START_REF] Maleki | Novel Leaking-Free, Green, Double Core/Shell, Palladium-Loaded Magnetic Heterogeneous Nanocatalyst for Selective Aerobic Oxidation[END_REF]. This catalyst can be used for promoting various reactions, including hydrogenation reactions [START_REF] Dehghani | Study of the effect of the ligand structure on the catalytic activity of Pd@ ligand decorated halloysite: Combination of experimental and computational studies[END_REF].

We are interested in using natural compounds for the catalytic purposes [START_REF] Sadjadi | Halloysite-based hybrids/composites in catalysis[END_REF][START_REF] Sadjadi | Halloysite nanoclay decorated with 2-amino pyrimidine functionalized poly glycidyl methacrylate: An efficient support for the immobilization of Pd nanoparticles[END_REF][START_REF] Sadjadi | Pd nanoparticles immobilized on the poly-dopamine decorated halloysite nanotubes hybridized with N-doped porous carbon monolayer: A versatile catalyst for promoting Pd catalyzed reactions[END_REF]. In this regard, we disclosed use of CS as a potential support for the stabilization of various catalytic species [START_REF] Sadjadi | Functionalized chitosan polymerized with cyclodextrin decorated ionic liquid: Metal free and biocompatible catalyst for chemical transformations[END_REF][START_REF] Esmaeilzadeh | Pd immobilized on hybrid of magnetic graphene quantum dots and cyclodextrin decorated chitosan: An efficient hydrogenation catalyst[END_REF].

On the other hand, the importance of industrial processes motivated us to design our catalysts for these types of processes. In this context, we proved the efficiency of natural clays for hydrofinishing of PAO under mild reaction condition [START_REF] Tabrizi | Efficient hydro-finishing of polyalfaolefin based lubricants under mild reaction condition using Pd on ligands decorated halloysite[END_REF][START_REF] Sadjadi | Combined experimental and computational study on the role of ionic liquid containing ligand in the catalytic performance of halloysite-based hydrogenation catalyst[END_REF]. In fact, in our previous works on the hydrogenation of PAOs, we applied halloysite nanoclay as a supporting material for the immobilization of Pd NPs. In the following of our research in this domain, we wish to report design and synthesis of a new CS-based catalyst for PAO hydrogenation. To approach high interaction of Pd nanoparticles (NPs) with support, CS, a ligand, selected from a series of appropriate ligand candidates by molecular simulation via Density-Functional Theory (DFT) method, was grafted on CS. Subsequently, the catalyst was synthesized, characterized and then applied for hydofinishing of PAO under relatively mild reaction condition. Compared to our previous reports [START_REF] Tabrizi | Efficient hydro-finishing of polyalfaolefin based lubricants under mild reaction condition using Pd on ligands decorated halloysite[END_REF][START_REF] Sadjadi | Combined experimental and computational study on the role of ionic liquid containing ligand in the catalytic performance of halloysite-based hydrogenation catalyst[END_REF], in which the ligand was synthesized on the support via multi-step procedures, in this work, the ligand was separately synthesized and then attached on CS. Moreover, the synthesis of the ligand was very facile and realized via a one-step procedure. Furthermore, in our previous reports we utilized clay as a support, while in this work we used CS that is more available and less expensive. In fact, in this new catalyst we aimed to benefit from the abundant heteroatoms of CS for improving stabilization of Pd NPs. Notably, in one of our previous works [START_REF] Karimi | Pd on nitrogen rich polymer-halloysite nanocomposite as an environmentally benign and sustainable catalyst for hydrogenation of polyalfaolefin based lubricants[END_REF], we prepared a composite of halloysite-polymer and then palladated it to achieve a catalyst. However, in this work we do not use composite and only make a chemical modification on CS via grafting of a ligand.

The justification for using CS is that CS is the most abundant biopolymer with outstanding features, such as low cost, biocompatibility, availability and abundant functionalities on its backbone that can be further functionalized.

Experimental

2.1.Reagents

The used reagents for the preparation of the catalyst were as follow: CS (Mw = 50000-190000, with deacetylation degree ≥ 75%, viscosity 20 cP for 1 wt.%, in 1% acetic acid), (3isocyanatopropyl) triethoxysilane (IS), 2-aminopyrimidine, Pd(OAc)2, toluene, NaBH4 and methanol (MeOH). All of the mentioned reagents were provided from Sigma-Aldrich and applied as received. For conducting hydrofinishing test and examining the activity of the as-prepared catalyst, PAO was first prepared by using 1-decene and AlCl3. These reagents were purchased from Merck Co. Germany and used without further purification.

2.2.Instruments and apparatus

PERKIN-ELMER-Spectrum 65 spectrometer, with KBr pellets was applied for conducting Fourier transform infrared (FT-IR) spectroscopy with scan time of 1 s and spectral resolution of 2 cm -1 .

To prepare the KBr pellets, 1 wt. % of the sample was used. Philips CM30300Kv instrument was used for recording Transmission electron microscopy (TEM) images of Pd/CS-IS-Py. To perform this analysis, the catalyst was deposited onto a carbon coated copper grid. Metal NPs size distribution has been estimated from the measurement of around 120 NPs found in arbitrarily chosen area of the images using Image J program. Energy dispersive spectroscopy (EDS) and mapping analysis were performed by MIRA 3 TESCAN-XMU. X-ray diffraction (XRD) pattern of the as-prepared catalyst was obtained over the range of 2θ = 5-80° via Siemens, D5000 diffract meter using Cu-Kα radiation as the X-ray source. To carry out Inductively Coupled Plasma (ICP-AES), Vistapro was applied.

To characterize the as-prepared PAO and measure the hydrofinishing yield, the following analyses were performed: 1 HNMR spectroscopy (Bruker DRX 400 MHz NMR spectrometer in deuterated chloroform at 25 °C), kinematic viscosimetry at 40 and 100 ℃ and viscosity index (KV 40 , KV 100 and VI were measured according to the ASTM D445 and ASTM D 97 test methods respectively), gel permeation chromatography (GPC Agilent 1100) using PS standards. To prepare IS-Py, IS (2 mL) and 2-aminopyrimidine (0.76 g, 0.8 mmol) were mixed under argon flow at 60 °C for 12 h. Then, the resultant ligand, IS-Py, was washed repeatedly with diethyl ether and dried in oven at 40 °C.

Decoration of CS with IS-Py ligand

To incorporate IS-Py ligand on CS, CS (2.1 g) was suspended in dry toluene (50 mL) under ultrasonic irradiation (power 150 W). Then, IS-Py (1.7 g) was added to the reaction media and the mixture was refluxed under Ar atmosphere at 110 °C overnight. At the end of the reaction, the precipitate was separated via centrifuge and washed with toluene several times. At the end, the synthesized support, CS-IS-Py, was dried in oven at 70 °C overnight.

Stabilization of Pd NPs on CS-IS-Py

To stabilize Pd NPs on CS-Py, Pd(OAc)2 (0.06 g) was dissolved in EtOH and gently dropped in stirring suspension of CS-IS-Py (2 g) in EtOH (60 mL). Stirring was continued for 120 min under Ar atmosphere. Then, NaBH4 solution (20 mL, 0.2 N in methanol) was prepared and introduced in a drop wise manner to the aforementioned mixture. After 2 h, the black precipitate was filtered and the obtained product, Pd/CS-IS-Py, was washed with methanol and dried in oven at 80 °C for 12 h, Figure 1. To determine Pd loading on Pd/CS-IS-Py, ICP analysis was conducted. The result confirmed low Pd content, ~ 1 wt.%, in the catalyst structure.

2.4.PAO synthesis

Synthesis of PAO through oligomelization of 1-decene was conducted according to our reported procedure [START_REF] Sadjadi | Rationalizing chain microstructure in the polyα-olefins synthesized by cationic AlCl3/H2O catalytic system[END_REF]. Prior to oligomerization, the reactor was purged with Ar at elevated temperature (75 °C) for 1 h. For oligomerization reaction, AlCl3 (5 g) was added to the reactor and then 1decene (500 mL) mixed with 0.3 mL deionized water (as electron donor, H2O/AlCl3 = 0.5 mol/mol) was fed and the reaction mixture was kept under stirring at 100 °C for 1 h. At the end of the reaction, the reactor was discharged and the obtained product was rinsed with a solution of 5 wt. % NaOH repeatedly. Then, the mixture was heated under -0.8 bar vacuum up to 250 °C to separate the obtained oil from the unreacted monomer, residual water and low molecular weight byproducts. Notably, the yield of the process was 86 %. According to the performed analysis (Figure 2), the synthesized PAO had Mn=1710 g/mol, PDI=1.4, KV 40 =245.5 cSt, KV 100 =32.1 cSt and VI=132.

2.5.Hydrofinishing of PAO

Hydrofinishing of PAO was carried out under relatively mild reaction condition. Briefly, PAO (10 g) and Pd/CS-IS-Py (0.5 g) were placed in a stainless steel reactor that was purged with dry nitrogen at 100 °C for 1 h. Subsequently, the reactor was sealed, the temperature was elevated to 130 °C and H2 gas with pressure of 6 bar was purged. The reaction mixture was stirred under the aforementioned condition for 6 h. Upon completion of the reaction, the reactor was cooled and discharged. Then, Pd/CS-IS-Py was separated via centrifuge and the yield of hydrofinishing was evaluated via 1 HNMR spectroscopy.

2.6.Computational details

All molecular simulations were undertaken using Gaussian 16 package [START_REF] Frisch | Gaussian 16[END_REF]. In this line, geometry optimizations were accomplished without any constraints via the BP86 functional of Becke and Perdew [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF][START_REF] Perdew | Density-functional approximation for the correlation energy of the inhomogeneous electron gas[END_REF]. To describe the electronic configuration of the system, the split-valence basis set, defined as Def2-SVP keyword in Gaussian, [START_REF] Weigend | Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy[END_REF] was employed for all studied atoms. The Pd atom was treated without pseudopotentials. To approve the nature of the stationary points as global minima, frequency calculations were also done. In all calculations, a 6-member ring with one dangling NH2 and two OH groups was utilized as support model. The designed ligand was lodged on this single ring.

The polarizable solvation model (SMD) variation of IEFPCM, developed by Truhlar et al. [START_REF] Marenich | Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions[END_REF], was used to define solvent effects in toluene, acetonitrile, ethanol, and tetrahydrofuran (THF) media, employing the same calculation level used in geometry optimizations.

Result and discussion

Simulation Results

reaction of some commercially available NH2 containing aromatic compounds with (3isocyanatopropyl) triethoxysilane (IS). In this regard, the type of aromatic cycle, either as heterocyle or not, the number of N atoms in the aromatic ring and the number of NH2 functionals in the compound were surveyed (see Figure 1). In the molecular simulations, the energy of Pd adsorption of the functionalized CS with designed ligands was computed and compared, since proper stabilization of catalyst moiety on the support structure is a fundamental issue in the evaluation of a heterogeneous catalyst [START_REF] Bahri-Laleh | Computational modeling of heterogeneous Ziegler-Natta catalysts for olefins polymerization[END_REF]. In fact, the higher adsorption energy is a sign of less catalyst leaching during targeted reaction. Adsorption energy of Pd stabilization on the designed supports containing 12 different amine compounds (Figure 1) is gathered in Table 1. It is clear that ΔEads and ΔGads have a similar trend. Notably, 6-member rings demonstrated higher adsorption energies (ΔEads <-37.5 kcal/mol) than 5-member ring counterparts (ΔEads > -36.3 kcal/mol), even with less N heteroatom in the cycle. Furthermore, we could not find a clear relationship between the number of N atoms in the amine structure and adsorption energy. Among the designed compounds, the combination of A7 amine, with a six member aromatic ring containing two N heteroatom and one NH2 group, with IS gave the best result with the highest ΔEads and ΔGads of -39.9 and -31.5 kcal/mol, respectively. Followed with this catalyst, the designed ones with A2, A6 and A10 amines beard the next highest ΔEads of -38.5 kcal/mol, which was close to it in A7 containing catalyst. Therefore, catalysts with A2, A6, A7 and A10 amines were selected as the next most proper ones, from energy point of view, and further investigation on the free volume calculation of the final catalysts was performed on these systems. To determine quantitatively the amount of free volume at the first sphere around catalyst center (Pd, where hydrogenation process takes place), SambVca2.1 package of Cavallo et al. [START_REF] Falivene | SambVca 2. A Web Tool for Analyzing Catalytic Pockets with Topographic Steric Maps[END_REF][START_REF] Falivene | Towards the online computer-aided design of catalytic pockets[END_REF] was employed to get steric map of the catalysts having high amount of ΔEads (<-38.0 kcal/mol). According to Figure 3 and Table 1, a free volume of 51.7 % was obtained around Pd centre for Pd/CS-IS-Py catalyst with A7 ligand, which was slightly higher than it in other catalysts. Indeed, this amount of Vfree renders good accessibility to the big PAO chains to approach catalyst active site and take part in the hydrogenation process. Concluding adsorption energy and free volume sections, catalyst with A7 amine was selected as the catalyst of choice, with the highest ΔEads and Vfree, which will impart a high efficiency to the designed catalyst. Next, to determine the effect of solvent type on the ΔEads and ΔGads of Pd(0) adsorption onto the selected ligand, the one containing A7 amine, structure optimizations were performed in toluene, acetonitrile, ethanol, and THF solvents. The obtained results of Table 2 emphasized the better performance of polar solvents, especially ethanol in stabilizing Pd metal on the designed catalyst.

Therefore, ethanol was selected as impregnation solvent in the catalyst synthesis section.

Catalyst characterization results

To verify incorporation of ligand and Pd NPs, EDS analysis was carried out. As shown in Figure FTIR spectroscopy was also exploited for confirming structure of Pd/CS-IS-Py. In this regard, FTIR spectra of CS and Pd/CS-IS-Py were recorded and compared, Figure 7. As known [START_REF] Esmaeilzadeh | Pd immobilized on hybrid of magnetic graphene quantum dots and cyclodextrin decorated chitosan: An efficient hydrogenation catalyst[END_REF], CS characteristic bands appear at 3396 cm -1 , 2918 cm -1 1114 cm -1 and 1645 cm -1 . Comparing the two spectra, it can be concluded that the absorbance bands of CS can also be observed in the FTIR spectrum of Pd/CS-IS-Py. Moreover, the new band, observed at 1730 cm -1 , can be assigned to the -C=O functionality in the structure of the ligand. It is worth noting that other characteristic bands of ligand, i.e., the representative bands of -Si-O and -C=N functionalities that are supposed to appear at 1114 cm -1 and 1645 cm -1 overlapped with the characteristic bands of CS. In Figure 9, the comparison of TG curves of CS and the catalyst are presented. As shown, CS TG curve exhibits two weight losses due to the loss of water (at temperature about 130 ºC) and degradation of CS at about 320 ºC. Regarding the TG curve of the catalyst, the loss of water and degradation of CS can also be observed. Moreover, the comparison of the two curves indicated slightly higher stability of the catalyst that can originate from the presence of silane-containing ligand (18 wt.%).

Hydrofinishing of PAO with Pd/CS-IS-Py

Performance of Pd/CS-IS-Py catalyst was appraised for hydrofinishing of PAO. The reason for selecting this chemical transformation is that it is an industrially attractive process, essential for improving the features of PAO. Mostly, harsh condition (high temperature and pressure) is required for hydrofinishing of PAO. Hence, we aimed to design a catalytic process for promoting this reaction under relatively mild condition. To this purpose, the main reaction variables, i.e.

hydrogen gas pressure, reaction temperature and catalyst dosage were first optimized, Table 3. To initiate the experiments, hydrofinishing has been performed under relatively low hydrogen pressure (6 bar) at 130 ºC with use of different Pd/CS-IS-Py loading (3-6 wt.%), Table 3. It was inferred that by increasing the dosage of the catalyst from 3 to 5 wt.%, higher yield of the desired product was achieved in lower reaction time. However, further increase of Pd/CS-IS-Py loading had no effect of the yield of the reaction. Hence, 5 wt.% Pd/CS-IS-Py was considered as the optimum catalyst amount. Motivated by this promising result, the effect of reaction temperature

22

was also appraised by performing the reaction at different temperatures (110-150 °C), Table 3.

The results indicated that reaction temperature could also affect the reaction yield and time.

Hydrofinishing at 110 °C furnished the hydrogenated product in lower yields, while increase of this parameter to 130 °C improved the yield of the reaction. Notably, further increase of the reaction temperature was not beneficiary for the process. Therefore, 130 ºC was selected as the optimum reaction temperature. Finally, it was investigated whether use of lower hydrogen pressure could promote hydrofinishing efficiently. Unfortunately, the results established that decreasing of the hydrogen pressure led to the decrease of the hydrofinishing yield. Therefore, the optimum hydrogen pressure was selected as 6 bar. Notably, hydrofinishing of PAO under optimum condition (T= 130 °C, t =6 h and P =6 bar) led to the hydrogenated PAO in 98% yield. TON and TOF for this reaction were 195 and 32 (h -1 ) respectively. The yield of the reaction was measured using 1 HNMR. As seen in Figure 10, in the 1 HNMR spectrum of the hydrogenated PAO, no vinylic and allylic peaks was observed, implying the hydrogenation of unsaturated bonds. Moreover, the color of the PAO changed from yellow to white. To check whether introduction of ligand on CS can improve its catalytic activity, a control catalyst, Pd/CS was prepared through immobilization of Pd nanoparticles on CS via the exact procedure used for the palladation of the catalyst. Afterwards, the catalytic performance of Pd/CS for the hydrogenation of PAO oil under the optimum reaction condition was appraised and compared with that of Pd/CS-IS-Py. According to the results, Pd loading of Pd/CS (0.08 wt.%) was lower than that of the catalyst (1 wt.%). Furthermore, the catalytic activity of Pd/CS was also lower than that of Pd/CS-IS-Py. In fact, using Pd/CS, hydrogenated PAO was achieved in 83%. The lower catalytic activity of this control catalyst can be attributed to its lower Pd loading, confirming the contribution of IS-Py ligand for anchoring Pd nanoparticles.

3.3.Reaction mechanism

According to the previous reports [START_REF] Tabrizi | Efficient hydro-finishing of polyalfaolefin based lubricants under mild reaction condition using Pd on ligands decorated halloysite[END_REF], the plausible mechanism for PAO hydrogenation can be depicted as Figure 11. In this process, -C=C double bond is first complexed to the Pd (0) center to form a π-complex. Notably, computational studies, reported by our group [START_REF] Tabrizi | Efficient hydro-finishing of polyalfaolefin based lubricants under mild reaction condition using Pd on ligands decorated halloysite[END_REF], suggested that in this step H-H bond is not cleaved. In the next step, hydrogen transfer occurs to form ethyl ligand on Pd. Finally, second hydrogen transfer from the metallic center happens to furnish the saturated PAO.

3.4.Comparative study

To appraise whether the performance of Pd/CS-IS-Py for PAO hydrogenation is comparable with other catalysts, its activity and recyclability were compared with some other catalysts that have been used for hydrogenation of 1-decene-derived PAO, Table 4. As shown, the activity and recyclability of the tabulated catalysts are almost similar. However, the merit of Pd/CS-IS-Py is its facile synthetic procedure and use of abundant and low-cost starting materials for the preparation of the support. Although halloysite nanoclay is also a natural compound, it is significantly more expensive than CS. On the other hand, as mentioned in the Introduction section, the synthetic protocols used for the preparation of the listed catalysts in Table 4 were multi-step and more complicated. aggregation was observed upon repeated use of the catalyst. In fact, the average particle size of Pd NPs in the reused catalyst (after five runs) was 5.5±0.3 nm, which is slightly larger than that of the fresh catalyst (5.4±0.1 nm).

To evaluate the stability of Pd/CS-IS-Py upon recycling, the FTIR spectrum of the reused catalyst after fifth run was compared with that of the fresh one, Figure 13. As shown, the FTIR spectrum of the reused Pd/CS-IS-Py is similar to that of fresh catalyst and showed all of the absorbance bands, approving the stability of Pd/CS-IS-Py.

Fresh catalyst

Reused catalyst

Hot filtration test was conducted. In detail, the hydrofinishing reaction was stopped after 1 h, Pd/CS-IS-Py was separated via centrifuge and yield of the reaction was estimated via 1 HNMR.

Then, that reaction was continued under catalyst-free condition. After 5 h, the reaction was halted and the yield of hydrofinishing was measured. Comparison of the yield of the reactions after 1 and 5 h established that in the absence of Pd/CS-IS-Py no reaction progress was observed. This result indicated that Pd NPs did not leach in the course of hydrofinishing and catalysis was heterogeneous.

Conclusions

In attempt to increase the utility of CS in industrial processes, CS was applied as an efficient and biocompatible support for the immobilization of Pd NPs and devising a catalyst for hydrofinishing of PAO. Prior to stabilization of Pd NPs, CS was decorated with a proper ligand, IS-Py, that has been prepared from reaction of IS and 2-aminopyrimidine. In this regard, computational study was used for selecting the most efficient amine compound among twelve potential candidates. Then, the ligand-decorated CS was palladated and utilized as a heterogeneous catalyst for hydrofinishing of PAO under relatively mild reaction condition. The results implied that under hydrogen pressure of 6 bar at 130 ºC, hydrogenated PAO has been achieved in 98% yield. Gratifyingly, the present catalyst could be easily recovered and reused for at least five consecutive runs with only slight loss of its catalytic activity and insignificant Pd leaching. Furthermore, hot filtration test approved heterogeneous nature of catalysis.

2. 3 .

 3 Catalyst synthesis: Pd/CS-IS-Py 2.3.1. Synthesis of IS-Py ligand

Figure 1 .

 1 Figure 1. The synthetic procedure of the catalyst and the investigated compounds for the synthesis

Figure 2 .

 2 Figure 2. GPC curve of the synthesized PAO

Figure 3 .

 3 Figure 3. Steric map of the xy plane and optimized structure of the designed catalysts with a) A2,

4 ABFigure 4 .

 44 Figure 4. A: EDS and B: elemental mapping analyses of Pd/CS-IS-Py

Figure 5 .

 5 Figure 5. A: SEM images of A: CS and B: Pd/CS-IS-Py

Figure 6 .

 6 Figure 6. A: TEM image of Pd/CS-IS-Py and B: particle size distribution curve of Pd NPs.

Figure 7 .Figure 8 .

 78 Figure 7. FTIR spectra of CS and the as-prepared catalyst.

Figure 9 .

 9 Figure 9. TG curves of CS and the catalyst.

Figure 10 .

 10 Figure 10. 1 HNMR spectra of A) crude PAO and B) hydrogenated PAO using reaction conditions

Figure 11 .

 11 Figure 11. The plausible mechanism for PAO hydrogenation in the presence of Pd/CS-IS-Py.

1: HNT decorated with ionic liquid 2 :Figure 12 .

 212 Figure 12. A: Investigation of the recyclability of Pd/CS-IS-Py for the hydrogenation of PAO

Figure 13 .

 13 Figure 13. Comparison of FTIR spectra of fresh and recycled Pd/CS-IS-Py

Table 1 .

 1 Vfree and binding energy, in terms of ΔEads and ΔGads, of Pd adsorption onto the functionalized CS with different amine compounds of A1 to A12, at BP86/def2-SVP level of theory. Reported energies unit is kcal/mol.

	Amine compound	ΔEads 1	ΔGads 2	Vfree (%)
	A1	-30.7	-24.5	-
	A2	-38.5	-29.9	49.5
	A3	-36.2	-26.7	-
	A4	-35.1	-27.2	-
	A5	-36.2	-26.7	-
	A6	-38.5	-28.9	50.0
	A7	-39.9	-31.5	51.7
	A8	-37.9	-28.1	-
	A9	-37.8	-28.0	-
	A10	-38.5	-29.9	51.2
	A11	-38.1	-29.5	49.6
	A12	-37.5	-27.8	-

1

ΔEads =ECatalyst-EPd-ELigand,

2

ΔGads =GCatalyst-GPd-GLigand

Table 2 .

 2 Binding energy, in terms of ΔEads and ΔGads, of Pd adsorption onto the CS-IS-Py ligand in different solvents, at BP86/def2-SVP level of theory. Reported energies unit is kcal/mol.

		Ethanol	Acetonitrile	Toluene	THF
	ΔEads	-41.0	-40.6	-39.5	-40.4
	ΔGads	-33.5	-32.9	-31.3	-32.6

Table 3 .

 3 Optimization of reaction variables for hydrogenation of PAO under Pd/CS-IS-Py catalysis.

	Run	Catalyst	Temp.	Pressure	Time (h) Yield (%)
		(w/w %)	(°C)	(bar)
	1	3	130	6	8	86
	2	4	130	6	8	91
	3	5	130	6	6	98
	4	6	130	6	6	98
	5	5	150	6	6	98
	6	5	110	6	8	90
	7	5	130	5	8	80
	8	5	130	4	9	78

Table 4 .

 4 Comparison of hydrogenation of PAO (derived from 1-decene) in the presence of

	various catalysts.

Acknowledgements

This work has been supported by the Center for International Scientific Studies & Collaboration (CISSC), Ministry of Science Research and Technology. Naeimeh Bahri-Laleh thanks MOLNAC (www.molnac.unisa.it) for its computer facilities.