Jihane Aouni
email: aouni@sanofi.com

Jean Noel Bacro
email: jean-noel.bacro@umontpellier.fr

Gwladys Toulemonde
email: gwladys.toulemonde@umontpellier.fr

Pierre Colin
email: pierre.colin@sanofi.com

Loic Darchy
email: loic.darchy@sanofi.com

Bernard Sebastien
email: bernard.sebastien@sanofi.com

Design Optimization for dose-finding trials: A review

Keywords: Adaptive trials, Design optimization, Dose selection, Patient allocation, Utility functions 1

Dose selection is one of the most difficult and crucial decisions to make during drug development. As a consequence the dose-finding trial is a major milestone in the drug development plan and should be properly designed. This article will review the most recent methodologies for optimizing the design of dose-finding studies: all of them are based on the modeling of the dose-response curve, which is now the gold standard approach for analyzing dose-finding studies instead of the traditional ANOVA/multiple testing approach. We will address the optimization of both fixed and adaptive designs and briefly outline new methodologies currently under investigation, based on utility functions.

1 Introduction / Optimal design considerations Selection of a dose (or doses) to carry into confirmatory phase III trials is one of the most difficult decisions that need to be made during drug development and it is generally agreed that the high attrition rate in phase III studies may, in part, be due to inadequate dose selection. Dose-finding studies therefore play a major role in drug development and should be carefully designed and conducted.

The scope of this review will be limited to the optimization of the design of parallel group phase IIB dose-finding studies, leaving dose escalation dose-finding studies out of the scope. We will also address the problem of the number of subjects per dose level (both in terms of total sample size and/or distribution of the subjects within treatment arms) in the dose-finding trial.

In addition, we will not consider the traditional ANOVA approach that uses contrast statistics to compare dose levels (e.g., differences between active doses and placebo); we will, instead, focus on the model-based approach, that is based on the estimation of the dose-response profile.

Dose-response modeling

The more recent approaches in dose selection have clearly moved the methodology from the multiple testing methods to dose-response "modeling" based methods such as the MCP-Mod method (Multiple Comparison Procedure and Modeling), see [START_REF] Bretz | Combining Multiple Comparisons and Modeling Techniques in Dose-Response Studies[END_REF] and Aouni et al. (2019a). Modeling the dose-effect relationship then becomes a necessary, key task in drug development; for this purpose, various dose-response models, linear or nonlinear, have become standard and are well known: for instance, linear, logistic, E-max (simple, Sigmoid or standardized form), quadratic, exponential and linear in log-dose models , see [START_REF] Pinheiro | Model-based dose finding under model uncertainty using general parametric models[END_REF], [START_REF] Pinheiro | Analysis of Dose-Response Studies-Modeling Approaches[END_REF]. These models are presented in Table 1. Throughout this publication, we will use d to refer to the dose under consideration and θ for the vector of parameters of the dose-response function f (d, θ) that represents the mean response for patients receiving dose d.

The linear model is the simplest one. Its simplicity is an advantage but it suffers from severe limitations [START_REF] Vandenberg | Hormones and Endocrine-Disrupting Chemicals: Low-Dose Effects and Nonmonotonic Dose Responses[END_REF]). An example of a non-linear dose-response model is the E-max model: E 0 is the placebo effect, E max is the asymptotic (associated to an "infinite dose") upper bound of effect, as compared to placebo, and ED 50 is the dose giving half of E max as change from placebo. A more complex form of this model is its Sigmoid version, where an additional parameter is included, ensuring greater shape flexibility, the "Hill" exponent (the slope), h, reflecting the shape of the dose-effect curve.

Other monotonic (as a function of dose) dose-response models can be considered: logistic, exponential and linear in log-dose models.

The quadratic model has the ability to describe a non-monotonic dose-response relationship, in a convex/U-shape (when β 2 > 0), or in a concave/umbrella-shape (when

β 2 < 0). Name f (d, θ) Linear E 0 + δ × d Emax E 0 + E max × d/(ED 50 + d) Sigmoid Emax E 0 + E max × d h /(ED h 50 + d h ) Log linear E 0 + δ × log(d + offset) Exponential E 0 × exp(d/δ) Quadratic E 0 + β 1 d + β 2 d 2 Logistic E 0 + E max /(1 + exp[(ED 50 -d)/δ])
Table 1: Examples of dose-response models.

Generalities on Optimal designs

In its widest sense, a design is defined by the doses considered in the trial, the total sample size and the relative allocation fraction to doses. Even if the number of doses and dose spacings are briefly discussed, we will essentially focus on optimization in the relative allocation fraction to doses so that, throughout this publication, design will be essentially defined by a vector ξ = (w 1 , ..., w K ) t , representing the relative allocation fraction to doses.

Several optimality criteria can be defined leading to designs such as A, C, D, E, T and G-optimal designs studied by [START_REF] Rady | Relationships among several optimality criteria[END_REF], for instance. We will only briefly present D and C optimality criteria, which are, in practice, the most useful in the optimization of dose-finding study designs. It is important to note that the optimality criteria that we will discuss are 'local', in the sense that the designs are optimal for a given model, and for the assumed model parameters: some aspects related to model misspecification will be discussed in section 2.3.

D-optimal designs maximize the determinant of the Fisher information matrix related to the parameters defining the dose-response model. This corresponds to minimizing the volume of the confidence ellipsoid for the vector of parameters of the dose-response model, see [START_REF] Davison | Statistical Models[END_REF] for instance. This criterion has the advantage of considering all model parameters, instead of focusing on one particular dose. This property can be described with the help of the equivalence theorem of [START_REF] Kiefer | The equivalence of two extremum problems[END_REF]): an interpretation of this theorem states that, in the context of linear regression, a D-optimal design will minimize the maximum variability of the mean response estimation.

The objective of C-optimal design, in the context of dose-finding studies, is to minimize the variance of a specific function of the parameters, this function being most often defined by a target dose to estimate. In practice, the optimization of the design is conducted using specific algorithms. We can mention the Fedorov-Wynn algorithm [START_REF] Fedorov | Theory of optimal experiments[END_REF][START_REF] Wynn | Results in the construction of D-optimum experimental designs[END_REF]). Also the 'DoseFinding' R-package [START_REF] Bornkamp | DoseFinding: Planning and Analyzing Dose Finding Experiments[END_REF]) provides a function to compute D-optimal designs.

Dose-finding based on a Pharmacometric approach

When it can be implemented, a pharmacometric approach based on exposure response modeling, instead of a doseresponse modeling, can be very efficient. The exposure response modeling approach consists in relating, through a parametric model, the plasma concentration (or a summary of the kinetics of the plasma concentration such as the AUC or the Cmax) of the compound to a biomarker of interest or the clinical endpoint. This approach is described in detail by [START_REF] Holford | Dose Response: Pharmacokinetic-Pharmacodynamic Approach[END_REF]. Exposure response modeling can be more informative than dose-response modeling because individuals differ in their exposure even if they receive the same dose. Furthermore, when exposure response modeling is based on biomarker response, it can benefit from the mixed modeling framework for longitudinal analyses. Once estimated, exposure response models can be used to select the doses through simulations (see [START_REF] Holford | Dose Response: Pharmacokinetic-Pharmacodynamic Approach[END_REF] for instance).

The exposure response approach is particularly useful when the recruitment of patients is challenging and the minimization of the study sample size is of major importance: this is the case for pediatric studies for instance.

In the context of exposure response analyses the same optimality criteria as those used for the dose-response models can be defined in the framework of mixed models. As an example, the search for D-optimal designs is implemented in PFIM software [START_REF] Bazzoli | Design evaluation and optimisation in multiple response nonlinear mixed effect models: PFIM 3.0[END_REF]). Those optimal designs can be used to minimize the number of blood samples per patient while maintaining sufficient precision in the estimates of parameters of the pharmacokinetics models: this is particularly important for designing pediatric pharmacokinetics studies.

Design Optimization of fixed design

In this section we focus on randomized, double-blind, placebo-controlled, fixed dose, parallel group dose-finding trials. In particular we will not consider cross-over trials or dose titration designs as, despite their appealing statistical properties, they are more complicated to implement for phase II dose-finding trials with a long treatment duration (see [START_REF] Ting | General Considerations in Dose-Response Study Designs[END_REF], for a short discussion).

Choice of dose range

The first task is to identify the range of dose of interest, that we will note [d m , d M ]: the dose will be selected within this dose range as well as the doses chosen for the dose-finding trial.

When d m = 0, the dose-finding trial includes the placebo arm as control (it includes the case when all the doses, including the null dose, are administered on top of a background medication). When it is non ethical to include a placebo arm, then d m > 0. In both situations there could be an active control arm in addition to the experimental drug dose arms.

The choice of d M is a little bit less straightforward: the choice of the upper limit of the dose range is based on safety considerations and is related to the available information on the Maximal Tolerated Dose (MTD) as identified in phase I [START_REF] Ting | General Considerations in Dose-Response Study Designs[END_REF]).

Number of doses, Dose spacings

The optimization of these two key components of trial designs are dependent on several aspects: the optimality criterion chosen, the primary objective of the trial (i.e. estimate the MED or EDp dose) and the true underlying doseresponse.

Concerning the Number of doses and dose spacing, these two points are addressed by [START_REF] Bretz | Practical considerations for optimal designs in clinical dose finding studies[END_REF]. Based on previous results obtained by [START_REF] Dette | Optimal designs for dose-funding studies[END_REF], where it is shown that the C-optimal (for MED and EDp definitions)

design can be explicitly calculated for common dose-response models with two or three model parameters, the authors give explicit formula defining the optimal design (giving the required doses and the percentage of patients allocated to those doses), for the Emax model, considering the C-optimality criterion (both MED and EDp definitions) as well as the D-optimality criterion. Their main findings are as follows (we assume that d m = 0): for the linear and log-linear models the optimal design requires two dose levels whereas it requires two or three dose levels for the Emax and exponential models.

In fact, the results provided by [START_REF] Bretz | Practical considerations for optimal designs in clinical dose finding studies[END_REF] or [START_REF] Dette | Optimal designs for dose-funding studies[END_REF] potentially identify the optimal design for the common dose-response models, but in practice, a direct application of these results for designing a dose-finding study can be difficult: a potential problem is that the optimal design only requires one intermediate dose between the placebo and the maximal dose (this can be problematic as it would then be impossible to explore alternative dose response relationships). To address this issue, one could envisage a more robust approach: optimize the design with respect to each model and keep the set of doses selected for the various models; again such an approach could not be really applicable as some of the intermediate doses identified can be 'odd' and not really feasible (this is also discussed in [START_REF] Bretz | Practical considerations for optimal designs in clinical dose finding studies[END_REF]); however we will address below, in more detail, robustness and model misspecification issues.

As a consequence a more pragmatic and realistic approach consists in choosing a priori the doses to be tested,

d 0 (= 0), d 1 , d 2 , ..., d K (= d M )
, and optimizing the designs in patients allocation only, through weights w 0 , w 1 , w 2 , ..., w K .

Given the maximal dose, generally the MTD, the doses could be chosen based on pharmacokinetics in minimizing the overlap in exposure between consecutive doses. Another simple approach, although justified by none of the optimality criteria, could be to choose equal spacing between consecutive doses. Other rules are proposed and discussed by [START_REF] Ting | General Considerations in Dose-Response Study Designs[END_REF].

Patient allocation

As mentioned above, we will now focus on the problem of optimizing patient allocation to the a priori chosen doses of the design: 0, d 1 , ..., d K . Mathematically, this optimization will consist in maximizing an objective function defined over a K-dimensional simplex defined as the set of vectors, ξ = (w k ) k=0,..,K with positive components such that k w k = 1. The various methods that will be discussed here have potentially the ability to reduce the sample size of the dose-finding study by 40% or 45% [START_REF] Bretz | Practical considerations for optimal designs in clinical dose finding studies[END_REF]) assuming that the true dose-response model is known and the true model parameter was used to perform the design optimization computation.

Optimality Criteria

For all optimality criteria that will be defined below, optimization of the design involves an objective function, Φ(ξ),

representing the efficiency of the design, that depends on the total information matrix:

M θ (ξ) = K k=0 w k I(d k ; θ)
where I(d k , θ) is the Fisher information matrix for one subject receiving dose d k and with θ as the dose-response model parameter value.

For the D-optimality criterion the optimal design maximizes the design optimality criterion:

Φ D (ξ) = log det(M θ (ξ)).
As mentioned above, choosing the D-optimality criterion, will result in maximizing the overall precision of the parameter estimates, as assessed by the volume of the confidence ellipsoid for the model parameters.

If we are more interested in obtaining precise estimates of some target doses, such as the MED or EDp defined above, then the C-optimality criterion must be chosen. This can be formalized as follows: the target dose is a function, ψ(θ), of the model parameters (explicit definitions have been given above for the MED and EDp for the Emax model), therefore the optimal allocation is defined by the design vector, ξ, that minimizes the variance var(ψ( θ)). Based on the delta-method [START_REF] Davison | Statistical Models[END_REF]) the objective is to minimize

var(ψ( θ)) ∂ψ ∂θ t var( θ) ∂ψ ∂θ
which is approximately equivalent to minimizing

∂ψ ∂θ t (N • k w k I(d k ; θ)) -1 ∂ψ ∂θ = ∂ψ ∂θ t (N • M θ (ξ)) -1 ∂ψ ∂θ
where N is the total sample size.

Therefore the application of C-optimality criterion consists in maximizing, in ξ, the design optimality criterion (ignoring N that has no influence on the optimization):

Φ C (ξ) = ∂ψ ∂θ t M -1 θ (ξ) ∂ψ ∂θ -1
.

Another example of C-optimality is used by [START_REF] Miller | Optimal Designs for Estimating the Interesting Part of a Dose-Effect Curve[END_REF]. Here the aim is to minimize the variance var(ψ( θ, d))

where ψ is a function of the parameter that also depends on the dose: ψ(θ, d) = f (d, θ) -f (0, θ), the mean difference versus placebo of the dose d. As defined above, the optimality criterion is then:

Φ(ξ, d) = ∂ψ(θ, d) ∂θ t M -1 θ (ξ) ∂ψ(θ, d) ∂θ -1
.

This optimality criterion can have a particular interest in the case when the dose-finding study is also the Proof of Concept (PoC) study. In this case, one of the primary objectives would be to characterize the effect of the maximal dose (d M ) and decide whether its effect is clinically relevant (subsequent analyses/studies would then be conducted) or whether the drug development must be stopped. In [START_REF] Miller | Optimal Designs for Estimating the Interesting Part of a Dose-Effect Curve[END_REF] the authors propose to primarily consider this maximal dose and define

Φ 1 (ξ) = Φ(ξ; d M )
as optimality criterion of a design. This criterion has an appealing interpretation: if the criterion value of a particular design is increased by x% as compared to the balanced design, then the balanced design requires an increase in sample size of x% to provide the same precision in estimates of interest as the design under consideration.

In many situations, focusing on the highest dose only is not appropriate. In this situation, a criterion quantifying the quality of a design should consider a dose interval [d L , d U ]: this interval could simply be the whole dose range [0,

d M ]
or a sub-interval, depending on the context and the study specific objectives. In this case, the authors propose to define

Φ 2 (ξ) = dU dL Φ(ξ, x) -1 dx -1
and use it as optimality criterion of a design.

Robustness, Handling of Uncertainty

The various methods of design optimization described above can lead to substantial reduction in total sample size (as compared to the balanced design for instance) but the calculations depend on the knowledge of the true dose-response function but also on the true values of the model parameters. Even though some knowledge is available at the beginning of phase II, the risk of misspecification in the dose-response models and/or the model parameter values cannot be neglected.

The problem of model misspecification is addressed by [START_REF] Bretz | Practical considerations for optimal designs in clinical dose finding studies[END_REF]. This is of great importance as, according to [START_REF] Dette | Optimal designs for dose-funding studies[END_REF], optimal design specifications are extremely sensitive to model misspecification (whereas they show some robustness against parameter value misspecification). Indeed, authors report potentially severe loss of efficiency of the nominal optimal design in case of model misspecification. As an example, if the local MED-optimal design for the log-linear model is used for the logistic model, then its asymptotic variance is approximately 100 times larger than the value that would be obtained if the local MED-optimal design for logistic model was used, see The first approach is a maximin approach: for each candidate design ξ its relative efficiency with respect to the optimal design of the models is computed:

ef f j (ξ) = C j (ξ)/C j (ξ * j )
where C j is the optimality criterion adapted to model j. Then the worst value, min{ef f 1 (ξ), ..., ef f m (ξ)}, is computed. And the robust optimal design is the design

ξ * = arg max ξ min{ef f 1 (ξ), ..., ef f m (ξ)}
that maximizes this worst value. This optimization can be unconstrained but maybe a more pragmatic approach would be to conduct the optimization within a set of pre-defined easily feasible (from a logistic point of view) designs.

The second approach proposed by [START_REF] Bretz | Practical considerations for optimal designs in clinical dose finding studies[END_REF] has more a Bayesian or model averaging spirit. It consists in assigning model weights (α j ) to each model j, reflecting the belief in each of the models (but the natural choice α j = 1/m, ∀j, can be a pragmatic option). And then the robust optimal design is the one that maximizes the weighted sum m j=1 α j log(ef f j (ξ)).

The problem of uncertainty in parameter values is addressed in [START_REF] Miller | Optimal Designs for Estimating the Interesting Part of a Dose-Effect Curve[END_REF] using a Bayesian approach. The authors propose a priori typical dose-response profiles (called 'scenario') and assign parameter values representative of (or consistent with) these scenarios, and they also assign prior probabilities α j to each of these scenarios, representing beliefs in these scenarios. And then using any of the design optimality criteria Φ, as defined above, the method consists in maximizing

m j=1 α j Φ(ξ, θ j )/Φ(ξ B , θ j )
where ξ B is the standard balanced design. This approach is more intuitive and less computationally demanding than specifying a full prior distribution π for the parameter values and computing the Bayesian average:

Φ(ξ, θ)/Φ(ξ B , θ)dπ(θ).
As a conclusion, we think that in case of much uncertainty in the dose response relationship the primary goal is to consider designs robust against model misspecification. For this purpose both the maximin or weighted sum of the log of relative efficiencies of [START_REF] Bretz | Practical considerations for optimal designs in clinical dose finding studies[END_REF] could be used; the maximin approach is maybe more conservative (as it protects against a worst case scenario) and easier to explain to a clinical team whereas the weighted sum approach is more flexible as it allows, in case of an interim analysis, to update the weights (using a Bayesian approach for instance) and then optimize the design of the subsequent patient cohorts based on the interim analysis results.

Approach based on utility functions

In the methods described above, the objective of the dose-finding trials were not really to identify an "optimal" dose (in a sense that would need to be clarified) but to identify the lower dose that achieves some targeted level of efficacy: this is the case for the MED or EDp doses discussed in the preceding sections. In order to introduce notions of optimality for the doses, it is necessary to take into account safety considerations. A relevant way to do so is to consider utility functions and to assign utility values to the doses.

Utility functions are generally introduced and defined within the decision theory framework. Utility describes the preferences of the decision maker and classifies/orders decisions. A decision theory result suggests that, in an uncertain (random) environment, all decision rules can be compared and classified using the expectation of a certain function, called the utility function. Various levels of complexity can be introduced within the framework of utility functions: economical/financial aspects could be accounted for in addition to efficacy and safety considerations.

Utility functions of various kinds (including a safety / penalty term penalizing high doses, rewards, global costs and/or costs per patient) have already been proposed and discussed in the context of drug development: we briefly summarize below the various approaches in the literature.

In [START_REF] Patel | A mathematical model for maximizing the value of phase III drug development portfolios incorporating budget constraints and risk[END_REF], utility functions are defined at the portfolio level: this work is beyond the scope of our review as the aim is to optimize the sample sizes allocated to the studies to the various drugs in a given portfolio in order to maximize the expected net present value. Other authors aim to optimize a given drug development plan using utility functions, sometimes defined as the expected net present value (as in [START_REF] Antonijevic | Optimizing Drug Development Programs: Type 2 Diabetes Case Study[END_REF], [START_REF] Kirchner | Utility-based optimization of phase II/III programs[END_REF], [START_REF] Patel | Designing phase II trials based on program-level considerations: a case study for neuropathic pain[END_REF]). In these works, the development plans are optimized with respect to phase II and phase III sample sizes, as well as Go/No-Go Decisions. In [START_REF] Foo | Designs to balance cost and success rate for an early phase clinical study[END_REF], the aim is to optimize (in terms of number of patients, number of blood samples per patient and corresponding sampling schedule) an early phase pharmacokinetics study with the objective of minimizing study cost (study failure has a specific cost that is taken into account in the optimization). For early phase studies, we can mention [START_REF] Hee | Decision-theoretic designs for small trials and pilot studies: a review[END_REF], where decisions theoretic based designs are reviewed.

The approach in [START_REF] Dragalin | Adaptive designs for dose-finding based on efficacy-toxicity response[END_REF] considers efficacy and safety simultaneously, described by two correlated binary variables. Defining a design as a probability measure, ξ, defined on a set of dose X, the authors define a penalty function

ϕ(x, θ, C E , C T ) = p 10 (x, θ) -C E (1 -p .1 (x, θ)) -C T
where p 10 (x, θ) is the probability of being a responder for efficacy without adverse reactions at dose x (and for model parameter value θ), p .1 (x, θ) is the marginal probability of having an adverse reaction (whatever efficacy) at the dose x, and the C E , C T are positive control parameters to be chosen. This penalty function ϕ can be interpreted as the inverse of a utility function: it decreases as the efficacy increases and it increases as the safety deteriorates. The authors define an integrated penalty

ϕ(x, θ, C E , C T ) ξ(dx),
and then a penalized information criterion:

I(x, θ)ξ(dx) ϕ(x, θ, C E , C T ) ξ(dx)
where I(x, θ) is the information matrix associated to the two, efficacy and safety, binary variables. A D-optimal design is a design that would minimize

det I(x, θ)ξ(dx) ϕ(x, θ, C E , C T ) ξ(dx) -1 = ( ϕ(x, θ, C E , C T ) ξ(dx)) M det I(x, θ)ξ(dx) -1
where M is the number of parameters for the models. Because of the penalty term, the optimal design will put more weight on doses for which the penalty (combining both safety and efficacy) is not too large.

A related approach, but with a slightly different objective, can be described as follows. Starting from a utility function, which can be similar to the inverse penalty in [START_REF] Dragalin | Adaptive designs for dose-finding based on efficacy-toxicity response[END_REF]:

U (d; θ) = U E (d; θ) × U T (d; θ),
where the U E , U T are, respectively, efficacy and safety related utilities (U E increases when the dose increases while U T decreases when the dose increases). Then, instead of searching for an optimal design that would give more weight to doses with high utility values (like in [START_REF] Dragalin | Adaptive designs for dose-finding based on efficacy-toxicity response[END_REF]), one could search for designs that provide the most accurate estimates of the utility values for some selected set or range of doses. This would then lead to the performance of design optimization using a similar approach to the computation of C-optimal designs. Assuming the utility function is sufficiently smooth and concave, the optimal dose d * is defined by

∂U (d * ; θ) ∂d = 0.
Then with a maximum-likelihood estimate of θ, θ, the optimal dose is estimated by

∂U ( d * , θ) ∂d = 0.
By application of standard differential calculus results, we conclude that d is a smooth function of θ, and that

∂d * ∂θ = - ∂ 2 U ∂d 2 -1 ∂ 2 U ∂d∂θ .
Using this relationship, we can then apply the methods for the determination of C-optimal designs with function d * (θ)

as a function determining the target dose (function ψ(.) of section 2.3). The application of this methodology is currently being assessed by the main author as part of her research work. But, to be applicable, such an approach relies on several assumptions: the parametric dose-response model should be perfectly specified and asymptotic normality should be achieved, whereas it is well documented that convergence issues might arise for the most complex models, like the logistic or sigmoidal Emax models, for instance.

Adaptive Design Optimization

In this section we will focus on the problem of conducting one or several interim analyses with the specific objective of improving the design of the dose-finding study and making better decisions at the end (choice of dose for phase III).

In this framework, the interim analysis results in updating patient allocation to the active doses and placebo.

Two-stage adaptive designs

We first focus on the simplest case where there is only one interim analysis: the first cohort of patients is enrolled according to a standard design (a balanced design for example) and then, based on the interim results, the patient al-location for the second part is optimized using one of the methodologies mentioned above. In practice, the following methods can be used:

(i) consider only one single dose-response model (Emax for instance), use the interim data to estimate the parameters and optimize the second part of the trial accordingly.

(ii) consider some uncertainty in model/parameter values: consider various dose-response models, estimate each model individually and then apply the maximin or the weighted efficiency approach of [START_REF] Bretz | Practical considerations for optimal designs in clinical dose finding studies[END_REF] or use the optimal Bayesian design approach of [START_REF] Miller | Optimal Designs for Estimating the Interesting Part of a Dose-Effect Curve[END_REF] in updating the the prior beliefs for each scenario. This latter approach is detailed in [START_REF] Bornkamp | Response-adaptive dose-finding under model uncertainty[END_REF].

When the data accumulated during the first part of the trial are included in the final analysis also, as in [START_REF] Miller | Optimal Designs for Estimating the Interesting Part of a Dose-Effect Curve[END_REF], there are other constraints to consider when optimizing the weights in the first stage. Indeed, in this case for each dose arm j its weight for final analysis, w j , is necessarily ≥ w j,0 = n j,interim N total , where n j,interim is the number of patients in arm j at interim analysis and N total is the total sample size of the trial at the final analysis.

Therefore, when optimizing the efficiency of the design, the w j weights are constrained to be ≥ w j,0 and not only to be ≥ 0. Once the final weights are obtained, the number of patients to be enrolled for the second part of the trial is

n j,(2) = N total × w j -n j,interim .
Other operational constraints can be handled quite easily. In case, for instance, there is a need to have a predetermined sample size or allocation rate in the placebo group, say w 0 , then the optimization would be conducted in the other allocation rates, with a new constraint on the sum of weights:

K k=1 w k = 1 -w 0 .
According to [START_REF] Miller | Optimal Designs for Estimating the Interesting Part of a Dose-Effect Curve[END_REF], the practical interest of performing an interim analysis with the only goal of optimizing the trial is questionable.

The first issue is related to the final statistical analysis: it is not fully clear to what extent statistical inference is disturbed by the adaptation following the interim analysis. This point is addressed using simulations by [START_REF] Miller | Optimal Designs for Estimating the Interesting Part of a Dose-Effect Curve[END_REF]: their results suggest a slight modification of the statistical properties of the estimates by the interim analysis.

In addition, the same authors conclude that there is no real evidence that design optimization following the interim analysis significantly improves the efficiency of the design. This lack of improvement brought by the interim design optimization is also noticed by [START_REF] Maccallum | Accounting for parameter uncertainty in two-stage designs for Phase II dose-response studies[END_REF].

However, such a statement must be mitigated:

(i) it is important to bear in mind that design optimization is not the sole objective of interim analyses. There are in fact many reasons that interim analyses are worthwhile: notably ethical reasons (ineffective or unsafe doses are dropped), or, for example, the need to bring information to the project team or management.

(ii) there is evidence (see section 3.2) that multiple interim analyses improve the performance of the design.

(iii) in [START_REF] Miller | Optimal Designs for Estimating the Interesting Part of a Dose-Effect Curve[END_REF] the interim analysis is conducted when the third of the patients are enrolled: the results might have been different if the interim analysis had been conducted when half of the patients were enrolled for instance; in fact, the timing of the interim analysis is a key parameter that drives the performance of the adaptive design, and it must be properly tuned and even optimized. Actually, the choice of the best timing should be made using simulations:

the relative efficiency, with respect to a fixed design, of the adaptive design would be assessed for adaptation occurring at various timing such as, for instance, when 20%, 30%, 40%, 50%, 60%, 70% of the patients are enrolled, and the timing associated with the highest relative efficiency should be chosen.

In [START_REF] Antonijevic | Impact of Dose Selection Strategies Used in Phase II on the Probability of Success in Phase III[END_REF], using the utility function approach, the authors compare various scenarios defined by the number of doses in the phase II studies, the type of design (fixed or Bayesian adaptive dose allocation design) with the objective of optimizing the probability of success in phase III and the expected net value: the authors point out the importance of allocating a sufficient number of patients in phase II and conclude that the Bayesian adaptive designs generally outperform fixed dose allocation designs.

Designs with multiple interim analyses

The principle is similar to the case of a single adaptation: at each of the interim analyses, the vector of relative allocation fraction to doses is updated based on new parameter estimates or belief updates in the dose-response models for the Bayesian design case. Because of the operational complexities induced by the interim analyses, such adaptive designs are more appropriate when the interval between patient enrollment and assessment and analysis of the efficacy criterion is very short. These types of adaptive designs are assessed through simulations in [START_REF] Dragalin | A Simulation Study to Compare New Adaptive Dose-Ranging Designs[END_REF].

The authors consider various types of adaptive Bayesian or non-Bayesian designs, as well as various optimality criteria, and assess the efficiency of these designs with 0 (corresponding to a fixed design), 1, 2, 4 and 9 equally spaced interim analyses. Their findings confirm the improvement brought by design adaptation, especially when both the number of doses and the dose range are large. In particular, increasing the number of adaptations leads to a gain in precision when selecting target dose(s).

Conclusion

The randomized, double-blind, placebo controlled dose-finding study is a key milestone in drug development. Its aim is not only to provide strong evidence of drug effect, but also to give guidance for choosing the most relevant dose for the later phase of a drug development. Therefore it is of primary interest to optimize this trial and maximize the likelihood of success of the trial.

However, the optimal designs obtained by computation can be difficult to implement in practice. In addition, the dose-finding study does not only consist in the characterization of the primary endpoint dose-response relationship:

another important objective is to provide some information on the safety profile for all doses used in the dose-finding study. In this case, design optimization should take into account both efficacy and safety endpoints, like in [START_REF] Antonijevic | Impact of Dose Selection Strategies Used in Phase II on the Probability of Success in Phase III[END_REF]. Current research conducted by the main author is focussed on this topic within the framework of utility functions [START_REF] Aouni | Utility-Based Dose-Finding in Practice: Some Empirical Contributions and Recommendations[END_REF]).

Table 7

 7

	of Dette

[START_REF] Bretz | Practical considerations for optimal designs in clinical dose finding studies[END_REF]

propose several methods, applicable to any optimality criterion: D-optimality or C-optimality related to either MED or EDp. All of these methods require, as a first step, identifying or proposing m candidate models f 1 , ..., f m , with their corresponding parameter values θ 1 , ..., θ m , and then for each of them, identifying their optimal design ξ * j .

Acknowledgement We would like to thank all the contributors and reviewers for their helpful comments and suggestions that led to an improved paper.

Conflict of Interest

The authors have declared no conflict of interest.


G A Bocan

H Breiss

S Szilasi

A Momeni

E M Staicu Casagrande

M S Gravielle

E A Sánchez

H Khemliche

Anomalous KCl(001) surface corrugation from fast He diffraction at very grazing incidence

Keywords: numbers: 34.35.+a, 34.50.-s, 79.20.Rf, 79.60.Bm

We present theoretical and experimental evidence of an anomalous surface corrugation behavior in He-KCl(001) for incidence along 110 . When the He normal energy decreases below 100 meV, i.e. He-surface distances Z > 2 Å, the corrugation unexpectedly increases up to an impressive > ∼ 85%. This is not due to van der Waals interactions but to the combination of soft potential effects and the evolution of He-cation and He-anion interactions with Z. This feature, not previously analyzed on alkali-halide surfaces, may favor the alignement properties of weakly-interacting overlayers.

The design of promising new materials in the form of thin film stacks brings along the challenge of producing such heterostructures to the highest quality and thus requires an accurate knowledge of the relevant features that rule quasi-epitaxial growth. One such feature is the substrate surface corrugation, which has been shown to drive the overlayer-substrate relative orientation even for weakly-interacting interfaces, such as organic/inorganic or transition metal dichalcogenides systems [1,2].

Grazing incidence fast atom diffraction (GIFAD) [3,4] is a highly sensitive and non-damaging probe of the projectile-surface interaction potential. Particularly well-suited for a thorough exploration of the surface electronic corrugation, it has already provided valuable structural information on a variety of systems [5][6][7][8][9][10][11], and its grazing geometry makes it compatible with in situ monitoring of thin film growth [12,13].

GIFAD occurs when keV-energy projectiles impinge on a surface along a low-index crystallographic direction, under axial surface channeling conditions [14]. The GI-FAD pattern arises from the combination of: a) interchannel interference, originated from the periodic array of channels, giving the Bragg peaks and b) intrachannel interference, originated from the corrugation of the interaction potential within a given channel, giving the rainbow peak as well as the supernumerary rainbows [15]. The result on the detection plane is a sequence of Bragg maxima modulated by the underlying intrachannel interference [15,16].

In this work we employ GIFAD to investigate the corrugation of the KCl(001) surface when probed with low normal energy 4 He 0 projectiles along the 110 channel. We present experimental data, information provided by the potential energy surface (PES) and simulations of the scattering dynamics, all of which put forward a counterintuitive and notorious increase of the surface corrugation with decreasing normal energy E ⊥ in the low E ⊥ < ∼ 60 meV range, i.e. for He-surface distances Z > 2.2 Å. Although long-distance surface corrugation increases have been experimentally observed by GIFAD for He/LiF(001) (∼ 12%) [17] and He/KBr(001) (∼ 45%) [18], this is, to our knowledge, the first time such an overwhelmingly large ( > ∼ 85%) increase is reported. Furthermore, we prove that the ruling factors of this puzzling feature are the evolution of the He-anion and He-cation interactions with Z and the continuance of the projectile-surface interaction during the whole scattering process (i.e. soft potential effects, SPE).

The experiments were performed at room temperature in a GIFAD setup similar to that of Ref. [13]. The KCl(001) sample was cleaved in air and immediately transferred to the ultra-high-vacuum chamber. Once mounted on a five-axis manipulator it was prepared insitu by annealing at about 650 K to ensure its cleanness and crystallographic order (corroborated from the well defined GIFAD patterns obtained [19]). The 4 He 0 beam was formed by neutralizing a He + beam with primary energies in the 300-600 eV range. Following the neutralization cell, a set of apertures of 0.4 × 0.09 mm 2 placed 0.36 m apart were used to reduce the beam divergence to less than 1 mrad. The two-dimensional angular distributions of scattered projectiles were collected by placing a detector formed by a Microchannel plate, a phosphor-coated screen and a CCD camera in the forward direction at a distance of 1.27 m. In Figure 1 we depict a GIFAD pattern for 600 eV 4 He atoms with E ⊥ =21 meV. The simulated pattern, that shows good accord with the experiment, was obtained with the Surface Initial Value Representation (SIVR) [20,21], a wellestablished semi-quantum approach to describe the scattering dynamics [22,23]. Note that collimation [21], surface defects and thermal vibrations may contribute to the elongation of the Bragg spots along the polar angle, though the latter two are not included in the simulations. Such polar spread does not affect our GIFAD-derived results, based on the elastic contribution to the intensity (θ f = θ i ± 0.035 • for the central Bragg order) [19]. The experimental corrugation, presented in Fig. 2(a), was obtained by processing the GIFAD patterns to get the effective interaction potential [7,19,24]. Since low-E ⊥ scattering takes place far from the surface, it should be ruled by a mild potential landscape, yielding a decreasing surface corrugation with decreasing E ⊥ but, intriguingly, we find that the corrugation increases with decreasing E ⊥ and this increase gets particularly sharp below 40 meV.

In order to gain more insight into this unexpected behavior we built a high-precision potential energy surface (PES) for the He/KCl(001) system based on Density Functional Theory (DFT), using the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [25] (henceforth the PBE-PES, further details in Ref. [22]). When the PES is averaged along the channeling direction (x), the axial potential V 2D (Y, Z) that actually rules the scattering process [26] is obtained, where Y is the position across the channelling direction. In Fig. 2(b) we show equipotential curves of the axial PBE-PES, from which the intrinsic corrugation for a given E ⊥ is obtained from the corresponding amplitude ∆Z = Z max -Z min .

Energy regions for which atoms scatter relatively far from the surface (Z > 2 Å) may bring along nonnegligible van der Waals (vdW) interactions, left aside in the PBE-PES. Hence, we built two additional PESs [27,28], where we replaced the PBE exchange-correlation functional by the rVV10 [29] and the DF2b86r [30] ones (henceforth respectively labeled vdW1 and vdW2), which include vdW interactions in a self consistent fashion though with different approaches. The intrinsic corrugations obtained with the three PESs give very similar curves (see Fig. 2(a)), in good accord with the trend of the experimental values, though about 25% below them. The intrinsic corrugation increase is thus not due to vdW interactions but should rather be explained through a thorough analysis of the axial PBE-PES [31].

The equipotential profiles across the 110 channel (see inset in Fig. 2(b)) have maxima over the Cl -rows. At mid-channel, over the K + row, there is a central minimum for E ⊥ < ∼ 60 meV and a local maximum for E ⊥ > ∼ 60 meV. The latter creates two off-center minima which migrate away from mid-channel with increasing E ⊥ . The equipotential curves show that the He-Cl - interaction (as obtained from DFT calculations) is always more repulsive than the He-K + one. The evolution of these interactions upon changing E ⊥ , particularly the slow sinking of the central minimum due to the vanishing He-K + interaction, makes Z min vary slower than Z max , explaining the intrinsic corrugation increase.

The residual discrepancy between the intrinsic and the experimental corrugations remains nevertheless intriguing, leading us to address the rainbow angle [32], a bright peak in the intrachannel contribution at maximal angular deflection [6,14,15,32]. The rainbow angle shows a wellknown high sensitivity to the projectile-surface interaction. While the corrugation affects mainly the intensity of the specular region to which only trajectories reflecting in the flattest regions of the axial PES contribute, the rainbow probes the steepest regions of the axial PES. For a given E ⊥ , and within a hard-corrugated wall (HCW) approach [33], we determined the PES-rainbow from the steepest slope of the corresponding equipotential curve, while the experimental rainbow was obtained from the processing of the GIFAD pattern [7,19,24].

In Fig. 3(a) we show the experimental rainbow angle together with those obtained from the three PESs considered. The increase observed in the experimental rainbow angle for decreasing low E ⊥ is analogous to the one obtained for the corrugation. The PESs rainbows show a mild increase, associated to the sinking central peak in the potential profiles, but they cannot reproduce the experimental rainbow trend for E ⊥ < 100 meV. The magnitude of the disagreement, once again unrelated to vdW interactions, strikingly contrasts with the fairly good accord discussed for the corrugation. Extracting information directly from the E ⊥ equipotential curves takes into account the projectile-surface interaction only at the turning point (HCW approach), thus neglecting the continuance of the projectile-surface interaction during the rest of the scattering process (i.e. soft potential effects (SPE)). We hence performed simula-tions of the scattering dynamics for each PES considered, by means of the semi-quantum SIVR approach [20,21] (henceforth labeled PBE-SIVR, vdW1-SIVR and vdW2-SIVR). The simulations were performed considering coherent illumination of a single channel (no interchannel interference) so that the rainbow angle can be directly obtained from the diffracted distribution. Strikingly, the PBE-SIVR rainbow reproduces the experiments in the energy range considered with an almost quantitative agreement. Including vdW contributions results in no significant differences, with vdW1-SIVR and vdW2-SIVR running slightly above PBE-SIVR.

The notorious difference observed between the PES and the SIVR rainbow angles, combined with the very good SIVR-Experiment accord, reveal in a quantitative manner the importance of SPE in central aspects of GI-FAD. These dynamical effects are naturally included in both the SIVR and the experimental rainbows (obtained by processing the GIFAD patterns) while the PES ones assume unrealistic broken-line trajectories.

Soft potential effects for the PBE-SIVR case can be visualized in Fig. 3(b) where equipotential curves of the axial PES are shown in 1:1 scale. For 3 < ∼ Z < ∼ 5 Å this plot shows a noteworthy shallow (< 7 meV) attractive region, due to polarization and lower-order attractive interactions within the DFT PES, which markedly alters the shape and corrugation of the equipotential curves. For high enough E ⊥ this well is inconsequencial, but low-E ⊥ projectiles traversing this region are strongly affected. This is illustrated in Fig. 3(b) for E ⊥ = 15 meV. Each equipotential curve deviates the trajectory according to the slope at their contact point so that a steeper slope favors E z → E y energy conversion. SIVR trajectories deviate from the HCW one upon entering the well region and also during their way out. In the asymptotic condition the projectile-surface interaction fades out and trajectories become linear but not parallel to the HCW trajectory. In the light of this analysis, the difference between the intrinsic and the experimental corrugations could also be due to SPE. Hence, from our SIVR simulations we selected trajectories contributing to the specular intensity. These trajectories are reflected either around the maximum or the minimum of the equipotential, with near-normal incidence and scattering directions. We can compute their SIVR accumulated phase difference (Eq. (10) in Ref. [20]) and obtain an effective corrugation from ∆Φ SIV R = 2 k ⊥ (∆Z) SIV R (depicted in Fig. 2(a)). The agreement with experiments is remarkable and shows that the PBE-SIVR calculation captures the factors determining the corrugation increase. From Fig. 2(a), the differences between the SIVR and PES curves, together with the SIVR-Experiment accord, provide a quantitative measure of SPE on the effective experimental corrugation obtained by means of the GIFAD technique.

In summary, we report a counter-intuitive increase of the corrugation and the rainbow angle for decreasing low-E ⊥ , by means of the GIFAD technique, in 4 He-KCl(001) along the 110 direction. This increase reaches an overwhelmingly large > ∼ 85% and is independently obtained from experiments and simulations. Its origin, unrelated to vdW effects, lies in the different evolution of He-cation and He-anion interactions, combined with SPE. As the latter do not strongly affect the corrugation dependence on E ⊥ , the experimental curve preserves the intrinsic corrugation trend, meaning that GIFAD provides reliable information on this surface feature. In contrast, the rainbow increase is almost fully explained in terms of SPE, strongly enhanced by the shallow well, whose steepest regions markedly deviate low-E ⊥ rainbow trajectories.

The same two factors are likely responsible for the corrugation increases at low E ⊥ reported for He/KBr(001) [18] and He/LiF(001) [17] (and expected for other alkalihalide surfaces) along this same channel. Particularly, for He/LiF, the much less spectacular corrugation increase is probably only due to SPE as the He-Li interaction is far vanished in the energy range considered. Also, we expect SPE to rule a marked rainbow increase feature, given the reported presence of shallow wells [34].

Thin alkali-halide layers (mainly NaCl and LiF) are often used to improve charge injection in organic electronic devices. The enhanced corrugation of KCl(001) at typical physisorption distances may favor the alignment properties of weakly interacting overlayers [1,2], improving their crystalline quality and thus the device performance. The measured corrugation may also be enhanced for probes either more reactive or with larger polarizability such as H or Ne [35][36][37]. Further investigations should address He GIFAD on other alkali-halide surfaces, aiming to maximize the intrinsic corrugation increase.

FIG. 1 :

 1 FIG. 1: GIFAD patterns for 4 He→ 110 KCl(001) with incidence energy E = 600 eV and E ⊥ = 0.021 eV. ϕ f (θ f ): azimuthal(polar) angular position of the scattered atoms.

FIG. 2 :

 2 FIG. 2: (a) Surface corrugation for 4 He→ 110 KCl(001) obtained from experiments, the three PESs considered and the respective SIVR simulations. (b) PBE-PES equipotential curves for 0 ≤ E ⊥ ≤ 200 meV. Inset: Surface geometry.

FIG. 3 :

 3 FIG. 3: (a) Rainbow angle for 4 He→ 110 KCl(001). (b) PBE equipotential curves (Only energies (meV) not shown in Fig. 2(b) are labeled) and trajectories contributing to the rainbow intensity (Dashed black line: HCW, Full black line: SIVR).

The authors acknowledge funding from ANPCYT (PICT-2017-1201; PICT-2017-2945) CONICET (PIP 112 201301 00386 CO; PIP 112 201501 00274 CO) and U.N. Cuyo (06/C590).


Navid Navidi-Izad

Hossein Hashemi

Borja Cascales-Miñana

Sylvie Régnier

Charles H Wellman

Thomas Servais

Colonial palynomorphs from the Upper Ordovician of northeastern Iran: spore "thalli", coenobial Chlorophyceae (Hydrodictyaceae) or cyanobacteria?

Keywords: acritarchs, spores, green algae, land plants, cyanobacteria

This study documents 'colonial' palynomorphs from the Upper Ordovician Ghelli Formation of northeastern Iran. The aggregates of organic-walled microfossils come from the Katian Armoricochitina nigerica and Ancyrochitina merga chitinozoan Biozones of this formation.

They can be classified as acritarchs and/or cryptospores, but they cannot be attributed to a particular biological group. Some specimens resemble 'thalli' of putative spores, such as Grododowon orthogonalis Strother 2017. Other clusters suggest an affinity to green algal groups, in particular with colonial chlorophyceaen algae, most probably of Hydrodictyaceae. Some specimens also show affinities with cyanobacterial groups. There are no arguments to relate these colonial palynomorphs to primitive land plants, but they are possibly produced by ancient green algal lineages that may have been experimenting with a subaerial existence.

Introduction

Unraveling the origin of the earliest land plants (embryophytes) and identifying their first fossil occurrence continues to be challenging (e.g. [START_REF] Gray | Major Paleozoic land plant evolutionary bio-events[END_REF][START_REF] Graham | Origin of land plants[END_REF][START_REF] Kenrick | The origin and early evolution of plants on land[END_REF][START_REF] Delwiche | The evolutionary origin of a terrestrial flora[END_REF][START_REF] Edwards | The early evolution of land plants, from fossils to genomics: a commentary on Lang (1937) 'On the plant-remains from the Downtownian of Wales[END_REF][START_REF] Wellman | The terrestrial biota prior to the origin of land plants (embryophytes): a review of the evidence[END_REF][START_REF] Gerrienne | Plant evolution and terrestrialization during Palaeozoic times-The phylogenetic context[END_REF]. For a long time, it was considered that mid-Silurian Cooksonia represented the oldest unequivocal macrofossil evidence for land plants [START_REF] Edwards | A vascular conducting strand in the early land plant Cooksonia[END_REF]. However, there are some recent claims, based on a variety of evidence, for an earlier origin of land plants (e.g. Strother et al. 2017;[START_REF] Morris | The timescale of early land plant evolution[END_REF]. It is generally considered that land plants emerged from a charophycean lineage of aquatic multicellular green algae (e.g. [START_REF] Graham | Origin of land plants[END_REF][START_REF] Karol | The closest living relatives of land plants[END_REF]).

However, [START_REF] Harrison | Development and genetics in the evolution of land plant body plans[END_REF] recently suggested that they may have originated from amongst a crust-forming terrestrial microbiome of bacteria, cyanobacteria, algae, lichens and fungi. Several green algal lineages have been proposed as the closest living relative to land plants, such as the Charophyceae, the Coleochaetophyceae and the Zygnematophyceae (see discussion in [START_REF] Gerrienne | Plant evolution and terrestrialization during Palaeozoic times-The phylogenetic context[END_REF].

Molecular clock evidence indicates origins of the different land plant lineages much earlier than the Silurian, at least in the Cambrian and possibly in the Precambrian [START_REF] Qiu | The deepest divergences in land plants inferred from phylogenomic evidence[END_REF][START_REF] Clarke | Establishing a time-scale for plant evolution[END_REF][START_REF] Zhong | Two new fern chloroplasts and decelerated evolution linked to the long generation time in tree ferns[END_REF][START_REF] Morris | The timescale of early land plant evolution[END_REF]. To date, however, no definite evidence of land plants exists before the Ordovician [START_REF] Kenrick | A timeline for terrestrialization: consequences for the carbon cycle in the Palaeozoic[END_REF][START_REF] Wellman | Palaeophytogeography of Ordovician-Silurian land plants[END_REF]Servais et al. in press). The microfossil record identifies the earliest cryptospores in the Middle Ordovician [START_REF] Strother | New evidence for land plants from the lower Middle Ordovician of Saudi Arabia[END_REF][START_REF] Rubinstein | Early Middle Ordovician evidence for land plants in Argentina (eastern Gondwana)[END_REF], while the first macrofossil records the oldest sporangia and small plant fragments in the Late Ordovician [START_REF] Wellman | Fragments of the earliest land plants[END_REF][START_REF] Salamon | Putative Late Ordovician land plants[END_REF]. Similarly, organic geochemistry investigations on biomarker compositions indicate a transition from green algae to land plants during the Early Palaeozoic, with the oldest reported biomarkers related to bryophytes from the Middle Ordovician of the Canning Basin, Australia [START_REF] Spaak | Environmental conditions and microbial community structure during the Great Ordovician Biodiversification Event; a multi-disciplinary study from the Canning Basin, Western Australia[END_REF].

By far the most significant amount of data on Ordovician land plants is provided by microfossil evidence. However, there is an ongoing discussion on the precise biological affinities of the dispersed spores and spore-like microfossils. Were these produced by genuine embryophytes or by organisms on the transition between ancient green algal lineages and land plants? Some authors suggest that some Cambrian 'cryptospores' were probably the desiccation-resistant spores of cryptogams belonging to the charophyte-embryophyte lineage (e.g. [START_REF] Strother | A classification scheme for the cryptospores[END_REF][START_REF] Strother | Middle Cambrian cryptospores and the origin of land plants[END_REF][START_REF] Taylor | Ultrastructure of some Cambrian palynomorphs from the Bright Angel Shale, Arizona, USA[END_REF].

On the other hand, most authors agree that cryptospores from the early Middle Ordovician (Dapingian) of Argentina described by [START_REF] Rubinstein | Early Middle Ordovician evidence for land plants in Argentina (eastern Gondwana)[END_REF] are derived from the earliest genuine land plants (embryophytes). To date, there is a large number of publications reporting such cryptospores from different palaeocontinents. So, for instance, the microfossil record shows occurrences from most parts of Gondwana and its periphery (e.g. [START_REF] Gray | Caradocian land plant microfossils from Libya[END_REF][START_REF] Vavrdová | Some plant microfossils of possible terrestrial origin from the Ordovician of Central Bohemia[END_REF][START_REF] Richardson | Late Ordovician and Early Silurian cryptospores and miospores from northeast Libya[END_REF][START_REF] Wellman | Fragments of the earliest land plants[END_REF]Wellman et al. , 2015;;[START_REF] Rubinstein | Cryptospore assemblages from the Ordovician/Silurian boundary in the Puna region, north-west Argentina[END_REF][START_REF] Steemans | Ordovician and Silurian cryptospores and miospores from Southeastern Turkey[END_REF][START_REF] Steemans | Origin and radiation of the earliest vascular land plants[END_REF][START_REF] Hérissé | Environmental changes reflected by palynomorphs in the early Middle Ordovician Hanadir Member of the Qasim Formation, Saudi Arabia[END_REF][START_REF] Rubinstein | Early Middle Ordovician evidence for land plants in Argentina (eastern Gondwana)[END_REF][START_REF] Mahmoudi | The first cryptospore assemblages of Late Ordovician in Iran, Ghelli Formation, eastern Alborz[END_REF][START_REF] Spina | Latest Ordovician (Hirnantian) miospores from the NL-2 well, Algeria, North Africa, and their evolutionary significance[END_REF][START_REF] Strother | Cryptospores from the Hanadir Shale Member of the Qasim Formation, Ordovician (Darriwilian) of Saudi Arabia: taxonomy and systematics[END_REF][START_REF] Rubinstein | The palynological record across the Ordovician/Silurian boundary in the Cordillera Oriental, Central Andean Basin, northwestern Argentina[END_REF][START_REF] Ghavidel-Syooki | Cryptospore and trilete spore assemblages from the Late Ordovician (Katian-Hirnantian) Ghelli Formation, Alborz Mountain Range, Northern Iran: Palaeophytogeographic and palaeoclimatic implications[END_REF][START_REF] Vecoli | Middle Ordovician cryptospores from the Saq-Hanadir transitional beds in the QSIM-801 well, Saudi Arabia[END_REF], including South China (e.g. [START_REF] Wang | Latest Ordovician cryptospores from southern Xinjiang, China[END_REF]), but also from Laurentia (e.g. [START_REF] Vecoli | Palynology of the Ordovician Kanosh Shale at Fossil Mountain, Utah[END_REF], Baltica (e.g. [START_REF] Vecoli | Cryptospore assemblages from the Upper Ordovician (Katian-Hirnantian) strata of Anticosti Island, Quebec, Canada, and Estonia: palaeophytogeographic and palaeoclimatic implication[END_REF][START_REF] Badawy | Late Ordovician (Katian) spores in Sweden: oldest land plant remains from Baltica[END_REF], Avalonia [START_REF] Richardson | Late Ordovician and Early Silurian cryptospores and miospores from northeast Libya[END_REF][START_REF] Wellman | Cryptospores from the type area for the Caradoc Series (Ordovician) in southern Britain[END_REF] and Siberia (e.g. [START_REF] Raevskaya | Cryptospores from the Katian (Upper Ordovician) of the Tungus basin: The first evidence for early land plants from the Siberian paleocontinent[END_REF]. This indicates that a flora appeared in the Middle Ordovician and already displayed a world-wide distribution by the end of the Late Ordovician.

The oldest known trilete spores come from the Upper Ordovician (Katian) of Saudi Arabia [START_REF] Steemans | Origin and radiation of the earliest vascular land plants[END_REF], which have been interpreted as potentially an indication of the occurrence of earliest vascular plants (i.e. tracheophytes; Wellman et al. 2015). However, as these only occur in low numbers and some extant bryophytes may also produce trilete spores, the earliest occurrences of this type of spore in the Late Ordovician does not indicate unambiguously the presence of vascular plants. The presence of stem group land plants 'at a bryophytic grade of organization' (cryptophytes) since the Middle Ordovician is thus established, whereas the presence of the first vascular plants possibly dates back to the Upper Ordovician with such plants only becoming common in the Late Silurian.

The Upper Ordovician (Katian-Hirnantian) Ghelli Formation of northeastern Iran also includes a good representation of typical Late Ordovician cryptospores and trilete spores. These land plant derived spores are present within well-preserved palynomorph assemblages that include acritarchs and chitinozoans (Ghavidel-Syooki 2016). The first cryptospore assemblage from Iran was documented by [START_REF] Mahmoudi | The first cryptospore assemblages of Late Ordovician in Iran, Ghelli Formation, eastern Alborz[END_REF] from this formation in the Khosh Yeilagh area, eastern Alborz Ranges, Iran. More recently, Ghavidel-Syooki (2016) reported a wellpreserved and diverse cryptospore assemblage from the Ghelli Formation in the Pelmis Gorge section in the Alborz Mountain Ranges, northeastern Iran.

In addition to typical Upper Ordovician cryptospores, the Ghelli Formation of northeastern Iran also includes early trilete spores, as well as clusters of palynomorphs that resemble colonial organisms, of possible algal or plant origin. The present study focuses on these organisms. Some specimens can be compared with the recently published spore 'thalli' of [START_REF] Strother | Ordovician spore 'thalli' and the evolution of the plant sporophyte[END_REF] from the Middle Ordovician Kanosh Shale at Fossil Mountain, Utah. [START_REF] Vecoli | Palynology of the Ordovician Kanosh Shale at Fossil Mountain, Utah[END_REF] first published planar clusters of 'cryptospores' with thick and smooth walls from the Kanosh Shale, before Strother in [START_REF] Strother | Ordovician spore 'thalli' and the evolution of the plant sporophyte[END_REF] described the new genus Grododowon, that the authors considered to represent a post zygotic growth phase in the life cycle of an ancient charophytic alga that was evolving in response to subaerial conditions. The question arises if the Iranian microfossils thus also indicate the presence of charophytic algae in the Upper Ordovician of the Gondwanan margin. Other specimens from Iran resemble material attributed to coenobial chlorococcales, in particular hydrodyctyacean algae, while others show morphologies that may relate them to cyanobacterial groups. Here, we aim to report the different morphotypes of 'colonial' microfossils from the Upper Ordovician Ghelli Formation and to discuss their nomenclature and possibly biological affinities.

Geological setting

The studied section (37° 21 69.12 N and 56° 93 83.80 E) of the Upper Ordovician Ghelli Formation is located in the northeastern Iran, about 40 km southwest of Bojnurd and ca. 57 km northeast of Jajarm (Fig. 1). [START_REF] Afshar-Harb | The stratigraphy, tectonic and petroleum geology of Kopet-Dagh region[END_REF] introduced the type section of the Ghelli Formation approximately 2.5 km north of Ghelli. The Lower Palaeozoic successions in this area comprise four formations, in ascending order, the Mila, Lashkarak, Ghelli, and Niur formations (Fig. 2). The Mila Formation is mostly composed of limestones attributed to the Middle and Upper Cambrian [START_REF] St Cklin | New data on the Lower Palaeozoic and Precambrian of North Iran[END_REF]. It is partially exposed at the investigated section and is overlain by the Lashkarak Formation consisting of dark green-grey siltstone, shale and thin layers of limestone and sandstone; attributed to the Early Ordovician [START_REF] Ghavidel-Syooki | Palynostratigraphy and palaeogeography of the Cambro-Ordovician strata in southwest of Shahrud City (Kuh-e-Kharbash near Deh-Molla), Central Alborz Range, northern Iran[END_REF][START_REF] Pour | Early Ordovician (Tremadocian) brachiopods from the Eastern Alborz Mountains, Iran[END_REF][START_REF] Zadeh | First record of the Ordovician fauna in Mila-Kuh, eastern Alborz, northern Iran[END_REF] on the basis of trilobite, brachiopod and acritarch biostratigraphy. The overlying Ghelli Formation has a thickness of about 900 meters. [START_REF] Afshar-Harb | The stratigraphy, tectonic and petroleum geology of Kopet-Dagh region[END_REF] divided it into three members including, in ascending stratigraphic order, volcanic rocks, shale and sandstone, and syndeposited mélange (sensu Afshar-Harb, 1979) of shales, siltstones, and sandstones (Fig. 3). The first palynological studies by [START_REF] Ghavidel-Syooki | Palynostratigraphy and paleobiogeography of Lower Palaeozoic strata in the Ghelli area, northeastern Alborz Range of Iran (Kopet-Dagh region)[END_REF] and [START_REF] Ghavidel-Syooki | Biostratigraphy and Palaeogeography of Late Ordovician chitinozoans from the north-eastern Alborz Range, Iran[END_REF] indicated that the Ghelli Formation can be attributed to the Upper Ordovician. The overlying Niur Formation consists mainly of dark grey shale, siltstone, sandstone and limestone with abundant corals, brachiopods, crinoids, and palynomorphs providing evidence for an attribution to the lower Silurian [START_REF] Afshar-Harb | The stratigraphy, tectonic and petroleum geology of Kopet-Dagh region[END_REF][START_REF] Ghavidel-Syooki | Latest Ordovician-Early Silurian chitinozoans from the eastern Alborz Mountain Range, Kopet-Dagh region, Northeastern Iran: biostratigraphy and paleobiogeography[END_REF].

Material and methods

More than 160 surface samples were collected from the Ghelli Formation; of which 105 more palynologically promising ones were processed in the palynological laboratories of the Kharazmi University (Tehran, Iran) and the UMR 8198 CNRSEvo-Eco-Paleo research unit of the University of Lille (France). Following standard procedures, palynomorphs were extracted from fine-grained samples such as shale and siltstone, using standard palynological techniques. This involved cold hydrochloric (HCl), cold hydrofluoric (HF) and hot hydrochloric (HCl) acid maceration to remove carbonates, silicates, and fluorosilicates, respectively. The organic residues were sieved through 15 and 20 m nylon meshes, respectively, at the universities of Lille and Kharazmi. No oxidation or alkali treatments were undertaken. All processed samples contain palynomorphs such as acritarchs, chitinozoans, scolecodonts and cryptospores in varying abundances and states of preservation. The present study concerns the description of the colonial palynomorphs. Most of these were recovered from the upper part of the Ghelli Formation (Fig. 3), i.e., the strata that can be attributed to the A. nigerica and A. merga chitinozoan Biozones [START_REF] Ghavidel-Syooki | Latest Ordovician-Early Silurian chitinozoans from the eastern Alborz Mountain Range, Kopet-Dagh region, Northeastern Iran: biostratigraphy and paleobiogeography[END_REF][START_REF] Ghavidel-Syooki | Palynology and isotope geochronology of the Upper Ordovician-Silurian successions (Ghelli and Soltan Maidan Formations) in the Khoshyeilagh area, eastern Alborz Range, northern Iran; stratigraphic and palaeogeographic implications[END_REF]Navidi-Izad;personal observation). All rock samples, palynological slides and residues are stored in the palaeontological collections of the Faculty of Earth Sciences of Kharazmi University, Tehran, except for the residues suffixed by ".s" that deposited at the UMR 8198 CNRS of the University of Lille, France.

Observations and Comparison

The palynological investigations in the Ghelli Formation confirm the presence of rich palynomorph assemblages containing acritarchs, chitinozoans, scolecodonts, cryptospores and trilete spores, that have partly been described in detail in previous studies [START_REF] Ghavidel-Syooki | Palynostratigraphy and paleobiogeography of Lower Palaeozoic strata in the Ghelli area, northeastern Alborz Range of Iran (Kopet-Dagh region)[END_REF][START_REF] Ghavidel-Syooki | Biostratigraphy and Palaeogeography of Late Ordovician chitinozoans from the north-eastern Alborz Range, Iran[END_REF][START_REF] Ghavidel-Syooki | Latest Ordovician-Early Silurian chitinozoans from the eastern Alborz Mountain Range, Kopet-Dagh region, Northeastern Iran: biostratigraphy and paleobiogeography[END_REF][START_REF] Ghavidel-Syooki | Palynology and isotope geochronology of the Upper Ordovician-Silurian successions (Ghelli and Soltan Maidan Formations) in the Khoshyeilagh area, eastern Alborz Range, northern Iran; stratigraphic and palaeogeographic implications[END_REF][START_REF] Mahmoudi | The first cryptospore assemblages of Late Ordovician in Iran, Ghelli Formation, eastern Alborz[END_REF][START_REF] Ghavidel-Syooki | Cryptospore and trilete spore assemblages from the Late Ordovician (Katian-Hirnantian) Ghelli Formation, Alborz Mountain Range, Northern Iran: Palaeophytogeographic and palaeoclimatic implications[END_REF], 2017a, 2017b). The presence of typical marine palynomorphs indicates that all samples are of marine origin. The high proportion of cryptospores in the assemblages (Ghavidel_Syooki, 2016; current study) suggests a terrestrial input from plants that inhabited land that was close by.

The assemblages from the Ghelli Formation also contain other palynomorphs, including aggregates of mostly sphaeromorph forms, some of them forming colonies (depicted in Plate 2) Importantly, different types of colonial palynomorphs are present in the studied assemblages (Plates 1, 2), most of them resembling microfossils attributed to colonial Hydrodictyaceae or cyanobacteria in the previously published literature.

Interestingly, Ghavidel-Syooki (2016, pl. VIII, fig. 12) illustrated a single specimen from the Ghelli Formation (re-figured in Ghavidel-Syooki 2017b, pl. I, fig. 12) of a colonial palynomorph that he attributed to Muzivum graziniskii. Ghavidel-Syooki (2016, fig. 2) classified this enigmatic palynomorph as an acritarch and indicated its presence in all samples investigated from the Ghelli Formation. Muzivum graziniskii was originally described from the Devonian of Poland by [START_REF] Wood | New Devonian coenobial chlorococcales (Hydrodictyaceae) from the Holy Cross Mountains and Radom-Lublin region of Poland: their palaeoenvironmental and sequence stratigraphic implications[END_REF], who considered the taxon to belong to the hydrodictyacean algae. In the present study, the presence in the Ghelli Formation of palynomorphs similar to Muzivum graziniskii is confirmed. For instance, the specimen illustrated in Plate 1, Fig. 8, typically resembles the planar, uni-layered coenobial sheets with a straight or undulatory outline, with the individual cells being thin-walled, primarily square and rarely rectangular, as described in the original diagnosis of Muzivum graziniskii [START_REF] Wood | New Devonian coenobial chlorococcales (Hydrodictyaceae) from the Holy Cross Mountains and Radom-Lublin region of Poland: their palaeoenvironmental and sequence stratigraphic implications[END_REF]. Other specimens are also very similar, but it remains difficult to attribute them to the same taxon with certainty, and they are here tentatively attributed to the hydrodictyacean algae or to algal clusters (Plate 1 and 2). The presence of this type of colonial microfossil in the assemblages examined suggests a fresh-water input into the marine environment.

Other specimens recovered herein resemble the recently described taxon Grododowon orthogonalis Strother 2017. Strother et al. (2017) described the genus Grododowon as a cluster of thick-walled, tightly adherent, smooth-walled spore-like dyads forming planar sheets. [START_REF] Strother | Ordovician spore 'thalli' and the evolution of the plant sporophyte[END_REF] considered the taxon as 'cryptospores', i.e. spore-like microfossils sensu [START_REF] Strother | A classification scheme for the cryptospores[END_REF] and [START_REF] Strother | Spore-like microfossils from Middle Cambrian strata: expanding the meaning of the term cryptospore[END_REF]. The specimens from the Ghelli Formation include dyad-like 'cryptospores' with thick and smooth walls. The dyads are generally about 10 to 15 m long and 5 to 8 m wide and they are arranged in orthogonal clusters similar to the type material from the Middle Ordovician of the USA [START_REF] Strother | Ordovician spore 'thalli' and the evolution of the plant sporophyte[END_REF]). The specimen illustrated in Plate 1, Fig. 1 appears to approach the diagnosis of Grododowon orthogonalis Strother 2017.

It displays a cluster with four dyad 'cryptospores' in one direction of the sheet and eight dyad cryptospores in the other direction. Other specimens (Plate 1) also fit in the diagnosis of Grododowon orthogonalis, although for some specimens, the identification remains problematical. It is interesting to note that the size of the planar sheets in the material investigated here is limited, ranging only from 45 to 55 µm. Larger sheets of 'thalli,' as illustrated by [START_REF] Strother | Ordovician spore 'thalli' and the evolution of the plant sporophyte[END_REF] from the Middle Ordovician of North America, are not retrieved from the material examined. A more delicate palynological processing technique, such as that used for 'small carbonaceous fossils' (SCF's) (e.g. [START_REF] Harvey | Sophisticated particle-feeding in a large Early Cambrian crustacean[END_REF] would potentially result in the recovery of larger sheets.

Other clusters of spore like microfossils ('cryptospores') are present in the same samples. The individual cryptospores display the typical morphology of cryptospores described from the Upper Ordovician (including from the Ghelli Formation) and can be considered as being related to genuine land plants (e.g. Plate 1, fig. 6).

On the other hand, several specimens of clusters of more or less spherical palynomorphs, usually between 5 and 10 µm in diameter, are depicted in Plate 2. They neither resemble cryptospores nor the spore 'thalli' described by [START_REF] Strother | Ordovician spore 'thalli' and the evolution of the plant sporophyte[END_REF]. They are also clearly different from the normal, loosely arranged clusters of unornamented (sphaeromorph) acritarchs that are very common in the Palaeozoic, such as clusters of Synsphaeridium spp. or Symplassosphaeridium spp., that were usually attributed to prasinophycean algae (e.g. Le [START_REF] Hérissé | Middle Ordovician acritarchs and problematic organic-walles microfossils from the Sag-Hanadir transitional beds in the QSIM-801 well, Saudi Arabia[END_REF]). The clusters from the Ghelli Formation (Plate 2) are usually arranged in structured 'colonies', clearly attached to each other, in regular patterns, with either thick or thin walls. All these specimens are here tentatively attributed to algal colonies.

The presence of planar sheets and cryptospore clusters, in association with isolated cryptospores derived from primitive plants, thus indicates a variation of microfossils that might reflect the presence of different types of green algal groups.

Sphaeromorph clusters, colonies, coenobia and 'thalli' in the Ordovician: acritarchs or cryptospores?

Aggregates of several or numerous identical organic-walled microfossils have been attributed to simple clusters, colonies, coenobia or 'thalli' in the palynological literature. For example, [START_REF] Vavrdová | Coenobial acritarchs and other palynomorphsfrom the Arenig/Llanvirn boundary, Prague basin[END_REF] mentioned 'coenobial' acritarchs in her Ordovician assemblages from Bohemia, as did Wood & Turnau (1996) by attributing their microfossils to 'coenobial chlorococcales. ' Strother et al. (1996, fig. 3, l) used the term 'cluster of spore like cells' for the problematic organisms of the Middle Ordovician of the USA. [START_REF] Foster | Eomerismopedia maureeniae n.g. n.sp., a chroococcacean cyanobacterium from the lower Ordovician Coolibah Formation, Georgina Basin, Queensland, Australia[END_REF] applied the term 'colonial palynomorph' to their microfossils from the Ordovician of Australia, whereas Le Hérissé et al. (2017, pl. 7), for example, used the terms cluster, colony, and coenobia to describe the different aggregates of palynomorphs from the Middle Ordovician of Saudi Arabia. Finally, [START_REF] Strother | Ordovician spore 'thalli' and the evolution of the plant sporophyte[END_REF] used the term 'thalli' for the spore-like structures from the Middle Ordovician of Utah, USA.

Such aggregates have been largely reported in publications concerning Palaeozoic palynology, and their nomenclature remains inconsistent. Sometimes, such 'clusters' have been attributed to the acritarchs, i.e. to palynomorphs of unknown biological affinity. [START_REF] Evitt | A discussion and proposals concerning fossil dinoflagellates, hystrichospheres, and acritarchs[END_REF] defined the informal group of the acritarchs to include small organic-walled microfossils of unknown biological affinity with various symmetries, shapes and structures. [START_REF] Evitt | A discussion and proposals concerning fossil dinoflagellates, hystrichospheres, and acritarchs[END_REF] did not exclude colonial palynomorphs in his original definition. However, several authors specifically excluded colonial organisms from the acritarchs. For instance, [START_REF] Fensome | Acritarchs and fossil prasinophytes: and index to genera, species and intraspecific taxa[END_REF] in the compilation of their catalogue of all acritarch taxa, did not include the colonial microfossils or aggregates of unknown organic-walled microfossils in their listings of taxa, and excluded them from their definition of the acritarchs. Subsequently some other authors further restricted the definition of the acritarchs to marine, single-celled microphytoplankton (see discussion in [START_REF] Servais | Some considerations on acritarch classification[END_REF] for further details). This led to different concepts of the definition of the acritarchs.

According to the restricted definitions after [START_REF] Evitt | A discussion and proposals concerning fossil dinoflagellates, hystrichospheres, and acritarchs[END_REF] original diagnosis, the microfossils presented here would not be acritarchs, as they are colonial and probably not marine. However, by following the original definition of the acritarchs by [START_REF] Evitt | A discussion and proposals concerning fossil dinoflagellates, hystrichospheres, and acritarchs[END_REF], all the colonial palynomorphs, including the specimens described in this study, can be classified as acritarchs, because their precise affiliation to a biological group is not proven. Additionally, [START_REF] Evitt | A discussion and proposals concerning fossil dinoflagellates, hystrichospheres, and acritarchs[END_REF] never mentioned that acritarchs must be marine organisms. As a result, all the specimens illustrated and discussed here fit the original diagnosis of the acritarchs by [START_REF] Evitt | A discussion and proposals concerning fossil dinoflagellates, hystrichospheres, and acritarchs[END_REF].

On the other hand, some of the illustrated specimens could be attributed to the cryptospores (or to clusters of cryptospores), because they also fit in some of the diagnoses of this group, as some specimens have spore-like morphologies. The term cryptospore, like the term acritarch, has also different definitions and conceptions. Cryptospores are organic-walled microfossils that resemble spores ('sporomorphs'), but that do not display all the morphological criteria to allow them to be attributed to land plants. It is important to consider that [START_REF] Richardson | Miospores, correlation and age of some Scottish Lower Old Red Sandstone sediments from the Strathmore region (Fife and Angus)[END_REF] first defined the cryptospores as "non-marine sporomorphs (non-pollen grains) with no visible haptotypic features such as contact areas or tetrad marks. Single grains or monads, 'permanent' dyads and tetrads are included" (Richardson et al. 1984 p. 116). In this original definition, the cryptospores are thus considered as being of terrestrial origin ('non marine sporomorphs'). In addition, [START_REF] Richardson | Late Ordovician and Early Silurian cryptospores and miospores from northeast Libya[END_REF] extended this definition to include single spores, naturally separated from dyads, which possess a circular contact area. Subsequently, [START_REF] Strother | A classification scheme for the cryptospores[END_REF] defined the cryptospores as "non-marine sporomorphs (non-pollen grains) without those typical visible haptotypic features such as trilete marks or furrows which characterized tracheophyte spores and pollen grains. Single grains or monad, 'permanent' dyads and tetrads are included as are sporomorphs separated from polyads which may or may not preserve contact area". Both [START_REF] Richardson | Late Ordovician and Early Silurian cryptospores and miospores from northeast Libya[END_REF] and [START_REF] Strother | A classification scheme for the cryptospores[END_REF] consider thus the cryptospores as 'non-marine sporomorphs' with no clear biological affinity. However, [START_REF] Steemans | Miospores evolution from the Ordovician to the Silurian[END_REF] emended Richardson et al.'s (1984) original definition by considering all cryptospores as derived from land plants: "Alete miospores (non-pollen grains) produced by primitive embryophytes. Single grains or monad, 'permanent' dyads and tetrads, and sporomorphs from polyads which may or may not preserve contact area". [START_REF] Steemans | Miospores evolution from the Ordovician to the Silurian[END_REF] thus considered that the cryptospores are dispersed spores of the earliest land plants (embryophytic).

It appears evident that all of the specimens illustrated here do not fit the diagnosis of the cryptospores by [START_REF] Steemans | Miospores evolution from the Ordovician to the Silurian[END_REF]. Nevertheless, some of the specimens described in this study could be considered as cryptospores in the sense of [START_REF] Strother | A classification scheme for the cryptospores[END_REF] and [START_REF] Strother | Spore-like microfossils from Middle Cambrian strata: expanding the meaning of the term cryptospore[END_REF], i.e. without a biological affiliation, and by accepting that the 'spores' can be assembled or colonial. It is also clear that the evidence that these microfossils are derived from land plants, and that they are non-marine, cannot be provided.

In summary, the specimens presented here could thus be classified as 'acritarchs' and, to some extent, as 'cryptospores' following the definition of the group by [START_REF] Strother | A classification scheme for the cryptospores[END_REF]. They could be named aggregates, clusters, colonies, coenobia (a synonym of colonies). The terminology 'thalli' should be applied with care, because it may suggest a relation to thallophytes, an abandoned biological category to group algae, lichens and fungi, and thus implying a biological affinity.

Possible biological affinities

The most interesting question is not the classification in artificial groups used by palynologists, but that about the origin of these organisms: what is the biological affiliation of these coenobia and planar sheets from the Upper Ordovician of Iran? Do these microfossils represent spores and clusters of spores of various primitive land plants or of different green algal groups?

In palaeobotany, the term miospore is a general term representing all fossil plant spores smaller than 200 m, regardless of whether they are cryptospores, isospores, microspores, small megaspores, pre-pollen or pollen grains. However, not only land plants produce spores, but also many algal groups (for example zygnemataceaen algae produce spores) or fungal organisms (fungal spores), that all have been described from palynological assemblages. As long as the biological affinities of these 'spores' is unknown, they can be classified as acritarchs, that can be single celled or colonial, and either of marine or fresh-water origin.

Moreover, some of the microfossils that are spore-like, but that do not display all the morphological characteristics that allow them be classified as land plant miospores, can be named 'spores.' Such 'spores' (or for some authors 'cryptospores') could be affiliated to different algal groups (or even fungi), and not be directly related to embryophyte land plants.

Such spores correspond to the definition of the cryptospores sensu [START_REF] Strother | A classification scheme for the cryptospores[END_REF], i.e. by excluding a biological affinity to embryophytes (sensu [START_REF] Steemans | Miospores evolution from the Ordovician to the Silurian[END_REF]. In this way, Strother et al. (2017) considered that the taxon Grododowon orthogonalis Strother 2017 can be classified among the 'cryptospores.' Indeed, the authors considered that these organisms are not produced by land plants but by algae close to the subphylum Charophyta, representing a "zygotic growth phase in the life cycle of an ancient charophytic alga."

In terms of biological affinities, some of the colonial palynomorphs from the Palaeozoic, including those presented herein, have been attributed in the literature to different groups of green algae, including the chlorophyceaen (Chlorophyta) and zygnemaphyceaen (Charophyta) classes. For many years, several authors have considered featureless, simple sphaerical acritarchs and clusters of such simple acritarchs as 'prasinophycean phycomata'(e.g. [START_REF] Colbath | Review of biological affinities of Paleozoic acid-resistant, organic-walled eukaryotic algal microfossils (including "acritarchs")[END_REF][START_REF] Guy-Ohlson | Green and Blue-Green algae. 7B-Prasinophycean algae[END_REF][START_REF] Hérissé | Middle Ordovician acritarchs and problematic organic-walles microfossils from the Sag-Hanadir transitional beds in the QSIM-801 well, Saudi Arabia[END_REF].

Other simple, featureless palynomorphs arranged in regular colonies, with multiple numbers of regularly-arranged cells, have been compared, and sometimes attributed, to the Hydrodictyaceae, a family of the Chlorococcales, that belong to the Chlorophyceae (a class of the Chlorophyta). The colonial Hydrodictyaceae are considered to be exclusively fresh water (e.g. [START_REF] Colbath | Review of biological affinities of Paleozoic acid-resistant, organic-walled eukaryotic algal microfossils (including "acritarchs")[END_REF][START_REF] Batten | Green and blue-green algae[END_REF]. They have been illustrated in many palynological studies, including from the Early Palaeozoic. Another enigmatic palynomorph is the brackish or freshwater colonial Botryococcus, a genus of green algae of the order Chlorococcales that is present in the fossil record since the Precambrian (e.g. [START_REF] Colbath | Review of biological affinities of Paleozoic acid-resistant, organic-walled eukaryotic algal microfossils (including "acritarchs")[END_REF].

The specimens from the Upper Ordovician of Iran mostly resemble these different groups of green algae belonging to the chlorophyceans, and some of them can be compared to modern hydrodictyacean algae. The 'spore-like thalli' described by [START_REF] Strother | Ordovician spore 'thalli' and the evolution of the plant sporophyte[END_REF] might be related to the class of the Charophyceae (and to a charophytic algal lineage), but they possibly also may belong simply to a group green algae, for example of the class of the Chlorophyceae, and in particular of the Order Chlorococcales.

In addition, other simple, colonial palynomorphs (acritarchs) have been related to the cyanobacteria. The problematic colonial palynomorph Gloeocapsamorpha prisca was described from the Ordovician Baltic Shale Basin of Estonia, where it is so abundant that it forms an 'oil shale,' exploited for its hydrocarbon content. Gloeocapsamorpha is currently considered to be an equivalent of the modern Entophysalidaceae "cyanobacteria", which are mat-forming and stromatolite-forming microorganisms [START_REF] Foster | Gloeocapsomorpha prisca Zalessky, 1917: a new study: Part I: taxonomy, geochemistry, and paleoecology[END_REF]). Most interestingly, [START_REF] Foster | Eomerismopedia maureeniae n.g. n.sp., a chroococcacean cyanobacterium from the lower Ordovician Coolibah Formation, Georgina Basin, Queensland, Australia[END_REF] described other colonial microfossils from the Lower Ordovician of the Canning Basin, Australia, and attributed them to the new species Eomerismopedia maureeniae. This taxon was considered by the authors to be another cyanobacterium. Morphological comparisons with modern material, in particular with the genus Merismopedia Meyen 1839, allowed [START_REF] Foster | Eomerismopedia maureeniae n.g. n.sp., a chroococcacean cyanobacterium from the lower Ordovician Coolibah Formation, Georgina Basin, Queensland, Australia[END_REF] to consider that Eomerismopedia to belonging to the cyanobacterial family Chroococcaceae and in particular the order Chroococcales. Some of the material recovered from the Upper Ordovician of Iran morphologically resemble this family of cyanobacteria, indicating that some of the recovered colonial palynomorphs may well be of cyanobacterial and not green algal origin.

Conclusion

In this paper, some colonial palynomorphs from the Upper Ordovician Ghelli Formation (northeastern Iran) are presented. The palynomorphs come from the Armoricochitina nigerica and Ancyrochitina merga chitinozoan Biozones (Katian in age) of this formation, of which the fossil content indicates typical marine environments, with a terrestrial input of cryptospores washed in from nearly land. The colonial (or coenobial) organic-walled microfossils can be attributed to the acritarchs, and some of them also to the cryptospores, depending on the definition used. However, they cannot be clearly attributed to any biological group with certainty. Some specimens resemble the possible 'thalli' of spores, similar to Grododowon orthogonalis Strother 2017. Such specimens may indicate evidence for a transition between algae and plants. However, they are not unambiguously related to land plants, and have been related to charophytic algae by [START_REF] Strother | Ordovician spore 'thalli' and the evolution of the plant sporophyte[END_REF].

Several other elements suggest an affinity to green algal groups, in particular with colonial chlorophyceaen algae, most probably of Hydrodyctyaceae. Others, however, may represent cyanobacterial groups as previously discussed by [START_REF] Foster | Eomerismopedia maureeniae n.g. n.sp., a chroococcacean cyanobacterium from the lower Ordovician Coolibah Formation, Georgina Basin, Queensland, Australia[END_REF].

In summary, the assemblages from the Upper Ordovician of Iran fit in [START_REF] Gray | Major Paleozoic land plant evolutionary bio-events[END_REF] Eoembryophytic phase (Middle-Upper Ordovician) of earliest liverworts-like plants. They also correspond to Strother et al.'s (2010) thalloid bryophyte landscape and/or the Gerrienne et al.'s (2016) Proembryophytic interval that was proposed to characterize the long transition from green algae to liverworts since the late Precambrian.

Figure 1 :

 1 Figure 1: Geographical map of northeastern Iran showing the location of the investigated section.

Figure 2 :

 2 Figure 2: Simplified schematic lithostratigraphic sketch (not to scale) of the Lower Palaeozoic succession in the investigated area.

Figure 3 :

 3 Figure 3: Stratigraphical column of the Upper Ordovician at the Ghelli section and stratigraphical distribution of the examined material. Chitinozoan biostratigraphy after Ghavidel-Syooki & Vecoli (2007), Ghavidel-Syooki et al. (2011) and unpublished results from the authors.

Figure 1 21

 1 Figure 1

Figure 2 22 1 Figure 3

 213 Figure 2

2

 2

Acknowledgments 13 Navid Navidi-Izad acknowledges a) the palaeontological team of the Evo-Eco-Paleo unit (Lille University -CNRS) for making his research stay possible and b) financial support received from the Iranian Ministry of Sciences, Research, and Techonology (MSRT). Thomas Servais thanks Van Mildert College, University of Durham, for financing a Prowse research fellowship. Hamed Bouzari Nezhad is acknowledged for his invaluable assistance during field work in Iran. This is a contribution to IGCP n°653 'The onset of the Great Ordovician Biodiversification Event.'