Hoai-Thu Thai
email: hoai-thu.thai@inserm.fr

Christine Veyrat-Follet

Nicole Vivier

Catherine Dubruc

Gerard Sanderink

France Mentré

Emmanuelle Comets

A mechanism-based model for the population pharmacokinetics of free and bound aflibercept in healthy subjects

Keywords: Population pharmacokinetics, target-mediated drug disposition, aflibercept, VEGF, MONOLIX Running head: Population PK model of aflibercept

What this paper adds" statement

What is already known about this subject Anti-angiogenic drugs have been developed as an effective therapeutic strategy for inhibiting tumour growth. However, their pharmacokinetics (PK) and their ligand inhibition properties have not been well characterised. The binding to a circulating target, such as vascular endothelial growth factor (VEGF), makes the PK of these drugs more complex.

What this study adds

The underlying mechanism of disposition of aflibercept, where a saturable and high-affinity binding of aflibercept to VEGF was adequately characterised by the Michaelis-Menten approximation of a target-mediated drug distribution (TMDD) model. To our knowledge, it is the first published mechanism-based population pharmacokinetic model for an anti-VEGF drug.

Summary

Aim: Aflibercept (VEGF-Trap), a novel antiangiogenic agent that binds to VEGF, has been investigated for the treatment of cancer. The aim of this study was to develop a mechanism-based pharmacokinetic model for aflibercept to characterise its binding to VEGF and its pharmacokinetic properties in healthy subjects.

Methods: Data from two phase I clinical studies with aflibercept administered as a single intravenous infusion were included in the analysis. Free and bound aflibercept concentration-time data were analyzed using a nonlinear mixed-effects modelling approach with MONOLIX 3.1.

Results:

The best structural model involves two compartments for free aflibercept and one for bound aflibercept, with a Michaelis-Menten type binding of free aflibercept to VEGF from the peripheral compartment. The typical estimated clearances for free and bound aflibercept were 0.88 L/day and 0.14 L/day, respectively. The central volume of distribution of free aflibercept was 4.94 L. The maximum binding capacity was 0.99 mg/day and the concentration of aflibercept corresponding to half of maximum binding capacity was 2.91 µg/mL. Interindividual variability of model parameters was moderate, ranging from 13.6 % (V max ) to 49.8 % (Q).

Conclusion:

The present PK model for aflibercept adequately characterises the underlying mechanism of disposition of aflibercept and its nonlinear binding to VEGF.

Introduction

Angiogenesis, the development of new blood vessels from pre-existing vasculature, participates in a variety of physiological processes and disease states [START_REF] Baka | A review of the latest clinical compounds to inhibit VEGF in pathological angiogenesis[END_REF]. Its critical role in tumour development and progression was established 30 years ago [START_REF] Folkman | Isolation of a tumor factor responsible for angiogenesis[END_REF]. This discovery brought a new effective approach called anti-angiogenic therapy to cancer treatment, which consists in limiting blood supply to tumours by preventing angiogenesis. The development of new agents has attracted many researchers' interest in the pharmaceutical industry. To date, the best characterised and most highly validated anti-angiogenic approach involves targeting the vascular endothelial growth factor (VEGF) pathway [START_REF] Holash | VEGF-Trap: a VEGF blocker with potent antitumor effects[END_REF].

VEGF is the most potent pro-angiogenic growth factor, promoting the formation of blood vessels which is required for both normal and neoplastic tissue growth [START_REF] Baka | A review of the latest clinical compounds to inhibit VEGF in pathological angiogenesis[END_REF][START_REF] Moreira | Vascular endothelial growth factor (VEGF) inhibition-a critical review[END_REF]. VEGF binds to two highaffinity receptors (VEGFR-1 and VEGFR-2) on endothelial cells. This binding activates the intrinsic tyrosine kinase activity of their cytodomains, initiating intracellular signalling. VEGF is expressed in a large variety of malignant tumours, such as tumours of breast, brain, lung and gastrointestinal tract [START_REF] Salven | Serum vascular endothelial growth factor is often elevated in disseminated cancer[END_REF]. Blockade of VEGF pathway is therefore an effective therapeutic strategy for inhibiting tumour growth [START_REF] Moreira | Vascular endothelial growth factor (VEGF) inhibition-a critical review[END_REF][START_REF] Salven | Serum vascular endothelial growth factor is often elevated in disseminated cancer[END_REF][START_REF] Ferrara | Vascular endothelial growth factor: basic science and clinical progress[END_REF].

Aflibercept (also called VEGF-Trap, Regeneron Pharmaceutics/ Sanofi-aventis research) is a novel anti-angiogenic agent that binds to VEGF with a 1:1 ratio and prevents it from interacting with its receptors. A recombinant fusion protein consisting of the second Ig domain of VEGFR-1 and the third Ig domain of VEGFR-2 fused to the Fc portion of human immunoglobin IgG1, it has a higher affinity for VEGF-A (K d in vitro = 0.5 pM) than current anti-VEGF monoclonal antibodies [START_REF] Chu | Aflibercept (AVE0005): an alternative strategy for inhibiting tumour angiogenesis by vascular endothelial growth factors[END_REF][START_REF] Lockhart | Phase I study of intravenous vascular endothelial growth factor trap, aflibercept, in patients with advanced solid tumors[END_REF][START_REF] Sternberg | Systemic chemotherapy and new experimental approaches in the treatment of metastatic prostate cancer[END_REF]. Aflibercept also binds to VEGF-B and Placental Growth Factor (PlGF), which may be advantageous in some settings, such as malignant ascites where PlGF may mediate vascular permeability [START_REF] Sternberg | Systemic chemotherapy and new experimental approaches in the treatment of metastatic prostate cancer[END_REF].

Based on the mechanism of action, this drug undergoes a target-mediated drug disposition (TMDD), a term used to describe the phenomenon in which drug is bound with high affinity to its pharmacologic target such that this interaction is reflected in the pharmacokinetic properties of the drug. A general PK model for drugs exhibiting TMDD has been developed by Mager et al [START_REF] Mager | General pharmacokinetic model for drugs exhibiting target-mediated drug disposition[END_REF][START_REF] Mager | Target-mediated drug disposition and dynamics[END_REF]. This model describes the elimination pathway of drug plasma concentrations as the combination of first-order elimination from the central compartment and specific target binding clearance followed by internalization of drug-target complex. It, as well, characterises the turnover of the target. The full TMDD model is complex and generally overparameterised. The more information we have about free drug, bound drug and the target, the more TMDD model components and parameters can be adequately identified, although it is yet unclear which elements should be measured to estimate all the parameters in the full TMDD model. In order to overcome this problem, several simpler forms of TMDD model were proposed [START_REF] Mager | Quasi-equilibrium pharmacokinetic model for drugs exhibiting target-mediated drug disposition[END_REF][START_REF] Gibiansky | Approximations of the target-mediated drug disposition model and identifiability of model parameters[END_REF]. There are mainly three approximations: quasi equilibrium (QE), quasi steady state (QSS) and Michaelis-Menten (MM). The QE approximation is based on the assumption that the drug-target binding is much faster than all other system processes. If the rate of elimination of complex is not negligible, the QE approximation is replaced by the QSS approximation assuming that the drugtarget complex concentration changes more showly than the binding and internalization process.

The MM approximation describes the system when the target concentration is small relative to the free drug concentration and the dosing regimens result in the target being fully saturated [START_REF] Gibiansky | Approximations of the target-mediated drug disposition model and identifiability of model parameters[END_REF].

Pharmacokinetics (PK) of aflibercept were investigated in healthy subjects after single intravenous (i.v.) doses of 1 to 4 mg/kg and single subcutaneous (s.c.) dose of 2 mg/kg, in two phase I clinical studies as part of the drug's clinical development. Both free and bound aflibercept concentrations were assayed. The objective of this analysis was to develop a mechanism-based PK model for aflibercept in order to characterise its binding to VEGF and its pharmacokinetic properties in healthy subjects. The influence of covariate was not assessed in this study due to the limited number of individuals and their healthy status

Methods

Study design

The data for the population PK analysis were collected from two phase I, monocentric and randomised studies which were both carried out in healthy male subjects to assess the PK of aflibercept. The studies were approved by the independent ethics committees (Pharma-Ethics, South Africa and Ethik-Kommission der Landesärztekammer Baden-Württemberg, Germany).

They are performed according to recommendations of the 18th World Health Congress (Helsinki, 1964) and all applicable amendments. The volunteers gave their written informed consent after full explanation of all procedures involved in the studies.

Study 1 was a placebo-controlled, single-dose, sequential ascending dose study. Forty eight subjects were enrolled in this study and equally divided into four groups: one group receiving placebo and three groups receiving a single dose of 1, 2 or 4 mg/kg of aflibercept respectively administered as a 1-hour i.v. infusion.

Study 2 was an open-label, single-dose, crossover study. Two groups of 20 subjects were included in the study. The first group received a single i.v. dose of 2 mg/kg in the first period, with a 2 month follow-up period, followed by a single s.c. dose of 2 mg/kg in the second period.

The second group received the s.c. dose first, then the i.v. dose.

The data from subcutaneous administration in the cross-over study (study 2) was removed from the analysis because this route of administration was not pursued in the subsequent clinical development. Moreover, in this study, a carry-over effect was found and should have been taken into account in the modelling, but this required the modelling of the subcutaneous route. To avoid it and to work on homogenous data, only the i.v. infusion data in the first period were used in the population analysis.

Blood sampling schedules

In study 1, blood samples (4 mL) were taken at the following times: predose, 1 (end of infusion), 2, 4, 6, 8, 12, 24 hours post-start of administration on day 1, then on days 8, 15, 22, 29, 36 and 43 at the same morning time corresponding to 2 hours after the start of infusion on day 1.

In study 2, blood samples (4 mL) were taken at the following times: predose, 1, 2, 4, 6, 8 hours post-start of administration on day 1, then on days 2, 3, 5, 8, 15, 29, 43 of each period at the same morning time corresponding to 2 hours after the start of infusion on day 1.

Assay method

For both studies, free and bound aflibercept plasma concentrations were measured in all samples collected collected at Regeneron Pharmaceuticals, Inc using enzyme-linked immunosorbent assay (Elisa) method. Blood samples were collected in tubes (containing 1mL of citrate buffer, sodium citrate, and 4.2 mg of citric acid) and were centrifuged at 2000 g for 15 minutes at room temperature. Plasma was stored at -20°C until analyzed.

In the assay of free aflibercept, the human VEGF165 was initially adsorbed to the surface of a polystyrene solid support to capture the free aflibercept in subject samples. A mouse monoclonal antibody (reporter antibody), specific to an epitope on the VEGFR1 domain of aflibercept, was then bound to the immobilized complex and an enzyme-linkedantibody (peroxidase-conjugated Affinipure goat antimouse IgG Fc-γ) was bound to the immobilized mouse monoclonal complex. A luminolbased substrate specific for peroxidase was added to achieve a signal intensity that was directly proportional to the concentration of free aflibercept. The limit of quantification of free aflibercept was 15.6 ng/mL. The calibration curves ranges from 100 ng/mL to 1.56 ng/mL in two-fold serial dilution. The limit of quantification of free aflibercept was 15.6 ng/mL. The inter-day accuracy and precision ranged from 92.24% to 103.09% and 1.05% to 16.18%, respectively. The intra-day accuracy and precision ranged from 106.36% to 109.90% and 9.56% to 13.68%, respectively.

In the assay of bound aflibercept, the human VEGF165 was replaced by a non-blocking goat antihuman VEGF antibody for capturing the bound aflibercept in subject samples. The rest of the procedure was similar as the one used in the assay of free aflibercept. The limit of quantification for bound aflibercept is 43.9 ng/mL. The calibration curves ranges from 100 ng/mL to 8.78 ng/mL in a 1.5 serial dilution. The inter-day accuracy and precision ranged from 93.85% to 110.19 % and 0.38% to 16.17%, respectively. The intra-day accuracy and precision ranged from 115.70% to 123.34% and 1.09% to 1.58%, respectively.

Population PK analysis

The population PK analysis was performed using the MONOLIX program (version 3.1) implementing the SAEM algorithm. The model control files were written using MLXTRAN script. The early concentrations of bound aflibercept were often below the limit of quantification, thus the censored data of bound aflibercept, representing 32.5% of data, were taken into account and used for model development using the extended SAEM algorithm implemented in MONOLIX as an exact maximum likelihood estimation method [START_REF] Samson | Extension of the SAEM algorithm to left-censored data in nonlinear mixed-effects model: Application to HIV dynamics model[END_REF]. In this algorithm, the leftcensored data are simulated in a right-truncated Gaussian distribution, instead of being imputed by the LOQ value or half of LOQ value. The data of bound aflibercept contains also some observed concentrations which were reported with values below the LOQ reported for most of the data (see bottom left plot in figure 4).

The database included a total of 56 subjects, with 36 subjects receiving treatments from study 1 and 20 subjects receiving i.v. infusions at the first period from study 2. With respect to the law of mass action, the concentrations of bound aflibercept were converted into equivalent concentrations of free aflibercept by multiplying them with 0.717, the ratio of molecular weights between free and bound aflibercept. The units of free aflibercept and bound aflibercept concentrations were µg/mL and µg.eq/mL respectively.

PK structural model

The following strategy was used to develop the model. Free aflibercept concentration-time data were first modelled alone. Then, bound concentration-time data were included for simultaneous modelling.

The structural model for free aflibercept was developed by testing the following models: twocompartment or three-compartment model with first order and/or MM elimination.

In the next step, we developed a structural model including bound aflibercept. The TMDD model with association and dissociation rate constants (k on and k off ), reduced approximate TMDD models and other simpler models with linear binding constant were used to describe the joint evolution of the two entities.

Statistical model

Denoting f the function describing the model, the statistical model for observed concentration C ij of subjects i for sampling time t ij is:

C ij = f(θ i ,t ij ) + ε ij
where θ i is the vector of parameters of subject i and ε ij is the residual error.

The errors ε ij are assumed to be independent and normally distributed with a null mean and a heteroscedastic variance σ² ij , which was modelled using a combined additive and proportional model:

σ² ij = (σ a +σ p f(θ i ,t ij ) )²
where σ a and σ p are additive and proportional coefficients of residual error model respectively.

Two alternative residual error models, proportional model (σ a = 0) and additive model (σ p =0) were also evaluated for the residual variability.

The interindividual variability on all parameters was modelled with an exponential model, e.g.

for CL:

CL i = CL.exp (η i,CL )
where η i,CL denotes the random effect in subject i, CL i the individual clearance parameter and CL the typical value of the population. The use of an exponential model implies a log-normal distribution for the parameters. The η's (e.g. η i,CL ) are zero mean random variables with variance ω² (e.g. for CL, ω²CL). The ω²'s represent the variance of the random effects; the elements of the interindividual variance-covariance matrix, Ω, was modelled as diagonal, e.g. assuming no covariance between the η's.

The following strategy was used for model development. First, the structural PK models were developed with a combined residual error model and interindividual variability on all parameters.

Then, the residual error models were evaluated for the selected structural model. Finally, the interindividual variability on each parameter was tested for significance, and non-significant variability component were removed one by one starting with the smallest and least significant estimate.

The log likelihood (LL) was computed using importance sampling. The likelihood ratio test was used to discriminate between nested models through the difference in log likelihood (-2LL). A pvalue of 0.05 was considered statistically significant. For non-nested models, the model selection was based on the Bayesian information criterion BIC. The better model is the one with a smaller value of BIC [START_REF] Bertrand | Comparison of model-based tests and selection strategies to detect genetic polymorphisms influencing pharmacokinetic parameters[END_REF].

Model evaluation

Internal evaluation of the model was based on goodness-of-fit (GOF) plots, including plots of observations versus individual and population predictions and plots of normalised prediction distribution error (NPDE) [START_REF] Brendel | Metrics for external model evaluation with an application to the population pharmacokinetics of gliclazide[END_REF]. A visual predictive check (VPC) was used to assess model predictive performance, based on the simulation of 500 datasets.

Results

Data

The database for population PK analysis contained data from 56 healthy subjects involved in the two phase I clinical studies. A total of 1476 concentrations were used for model building: 732 concentrations of free aflibercept and 744 concentrations of bound aflibercept of which 242 (32.5 %) were below the quantification limit (LOQ= 43.9 ng/mL or 0.0314 µg.eq/mL).

The pooled concentrations of free aflibercept and bound aflibercept plotted versus time are presented in Figure 1.

The time-course of free aflibercept after i.v. infusion of doses of 1, 2 and 4 mg/kg in semi-log scale suggests a bi-exponential decline with a rapid phase of distribution followed by a prolonged terminal phase of elimination regardless of the dose. The time-course of bound aflibercept suggests a saturable binding phenomenon with the same observed plateau for two higher doses.

The peaks of complex occurred sooner for dose 1 mg/kg and later for doses 2 and 4 mg/kg (around 21 days, for both doses).

Free aflibercept modelling

Several models were proposed to characterise the kinetics of free aflibercept. In the twocompartment model with first-order elimination, the distribution of post-hoc clearance for 3 doses showed a dose-dependent clearance, confirming the nonlinear disposition already observed during non-compartmental analysis. The two-compartment model with MM elimination was then developed and showed a better description of data with a decrease of 33.61 point in -2LL value.

The three-compartment models with first-order or MM elimination were also evaluated. They provided a slight decrease in -2LL compared to the two-compartment models but the fitting did not show a significant improvement. The two-compartment model with MM elimination was therefore selected as an adequate structural model for free aflibercept.

The best residual error for this model was the combined additive and proportional error. The interindividual variability on K m was not significant.

Free and bound aflibercept modelling

The binding of aflibercept to VEGF and the nonlinear kinetics of free aflibercept suggests that this drug has target-mediated drug disposition properties. Therefore, the TMDD model with association and dissociation rate constants (k on and k off ), reduced approximate TMDD models and other simpler models with linear binding constant were tested. The MM approximation of TMDD model developed by Gibiansky and colleagues [START_REF] Gibiansky | Approximations of the target-mediated drug disposition model and identifiability of model parameters[END_REF] was found the best approach to describe the kinetics of both free and bound aflibercept, while others did not fit the concentrations of bound correctly and/or had overparameterization issues. In all of these models, the dissociation rate k off was very small compared to other constants and very badly estimated.

The nonlinear central and peripheral binding of this MM approximate TMDD model were tested and compared. The nonlinear peripheral binding model was found to better describe the data than the central binding one with a smaller value of BIC (2160 vs 2703) and therefore retained as the structural model for free and bound aflibercept.

The final structural model relating free and bound aflibercept is shown in Figure 2.

Free aflibercept in plasma distributes first to tissues then binds to VEGF to form the complex.

Binding to VEGF follows the law of mass action and can be characterised by a nonlinear form with MM constants (V max , K m ). The bound aflibercept (complex) is assumed to be directly eliminated through internalisation (k int ) and not through the dissociation rate constant (k off ) which gives back free aflibercept and free VEGF. The concentration of free aflibercept in central compartment (C p ), in tissue compartment (C t ) and the concentration of bound aflibercept (C b ) are described by means of differential equations:

(1) (2)
(3) in which k el is the first order elimination rate constant of free aflibercept (day -1 ) from central compartment, k tp and k pt are the first order rate constants between central and peripheral compartment (day -1 ), k int is the first order rate constant of bound aflibercept internalization (day -1 ), V max is the maximum binding capacity (mg/day), K m is the concentration of free aflibercept corresponding to half of maximum binding capacity (µg/mL), V p is the central volume of distribution of free aflibercept (L), V t is the peripheral volume of distribution of free aflibercept and V b is the volume of distribution of bound aflibercept (L). More details concerning the derivation of the model equations and its relationship to other TMDD models proposed in the literature can be found in the Appendix.

C k C K C V V dt dC C K C V V C k V V C k dt dC V V C k C k k Input dt dC . . 1 . 1 . . . . . ). ( int max max - + ⋅ = + ⋅ - - = + + - =
In this model, the volume of bound aflibercept (V b ) and the maximum binding capacity (V max ) were however highly correlated with a correlation coefficient higher than 0.95, suggesting identifiability issues. In order to prevent this problem, the value of V b was assumed to be equal to the value of central volume of distribution of free aflibercept (V p ), corresponding to the hypothesis that free and bound aflibercept are distributed within the same space in the plasma compartment. To do so, V b was fixed to the population value of V p with no interindividual variability, instead of being estimated.

The best residual error model of bound aflibercept was the proportional model, while that of free aflibercept was the combined additive and proportional error model found in the first step. A combined error model was also attempted; however, the contribution of the additive error was negligible and was therefore discarded. The interindividual variability on the internalization rate constant (k int ) was found to be small (11.9%) and badly estimated. The likelihood ratio test demonstrated that removing this variability from the model did not significantly change the fit.

The estimated parameter values for the best model of free and bound aflibercept are presented in Table 1.

All model parameters were well estimated with relative standard errors (RSE) < 20%. The population estimate of clearance for free aflibercept was 0.88 L/day and internalisation rate for bound aflibercept was 0.028 day -1 , resulting in a clearance of 0.14 L/day. The central volume of distribution for free aflibercept was 4.94 L. The maximum binding capacity was 0.99 mg/day and the estimated concentration of free aflibercept corresponding to half of maximum binding capacity was low (2.91 µg/mL).

Interindividual variability for random effects associated with model parameters was moderate, ranging from 13.6 % (V max ) to 49.8 % (Q). Residual errors were low for both free aflibercept and bound aflibercept with proportional parts of 17.1 % and 12.6 % respectively. The additive error of free aflibercept was also small compared to the observed concentrations. Figure 5 presents the visual predictive checks for both free and bound aflibercept. The 80% prediction interval and median were obtained by simulating 500 datasets under the final model; in addition, we obtained the 90% prediction intervals around the median and the upper and lower boundaries of the interval. Concentrations lower than the LOQ simulated for free aflibercept at the end of treatment were replaced by the LOQ (0.0156 µg/mL) in order to obtain a lower limit for the VPC plots in semi-logarithmic scale. The slight misfit for free aflibercept seen at the last time-point may be a consequence of the number of concentrations near the LOQ, since the median and upper boundary appear to be very well predicted by the model. For bound aflibercept, the model slightly underpredicted median concentration from day 21 onwards but correctly predicted the variability. Thus, the model described reasonably well the observed concentrations of both free and bound aflibercept.

Discussion

In this study, we report the development of a mechanism-based model to characterise the population PK of aflibercept after a single i.v. infusion of a 2 to 4 mg/kg dose to 56 healthy subjects. Free and bound aflibercept plasma concentrations were used in the modelling.

We first modelled free aflibercept concentrations alone. Using standard compartmental PK models, we observed a decrease of clearance for higher concentrations. This is consistent with previous results from noncompartemental analyses, which demonstrated that the volume of distribution of free aflibercept was low and its clearance was dose-dependent. When dealing with saturable kinetics, the most common model to consider is MM elimination kinetics, which was retained as the final model for free aflibercept during model development. Although such a model has been successfully applied to describe nonlinear PK of free aflibercept, it does not well represent the underlying molecular events such as target binding, internalization and degradation of the drug. However, the modelling of free aflibercept helped us get general information about the model structure: two compartment kinetics and dose-dependent clearance. More complex models combining linear and nonlinear elimination were also tested. They adequately fitted the data but were not retained as the final model of free aflibercept according to the statistical criteria.

The next step of modelling was then to add bound aflibercept data, which serves as a marker of efficacy by representing the amount of VEGF bound to aflibercept. The mechanism of action of this drug suggested the use of TMDD model. Among several types, the MM approximation of TMDD model developed by Gibiansky et al proved the most appropriate to reflect the kinetics of aflibercept in our study as free aflibercept concentrations were much higher than target concentrations and their binding to VEGF resulted in a fully saturated target. The relatively low level of free circulating VEGF was confirmed by the range of 19-47 pg/ml of plasma VEGF levels observed in healthy subjects in a meta-analysis [START_REF] Kut | Where is VEGF in the body? A meta-analysis of VEGF distribution in cancer[END_REF]. In order to model simultaneously both free and bound drug, the MM equation of bound aflbercept was added into the system of differential equations. Free aflibercept is therefore eliminated through two pathways: non saturable elimination from central compartment (k el ) and a specific and saturable binding to VEGF, followed by internalization (k int ) of bound aflibercept, which is its dominant elimination pathway. In addition, the formation of bound aflibercept could occur predominantly in central or peripheral compartment. The final PK model was an approximate TMDD model involving two compartments for free aflibercept and one compartment for bound aflibercept, with MM-type binding of free aflibercept to VEGF from the peripheral compartment.

The first-order dissociation rate constant (k off ) of the complex to give back free aflibercept was assumed to be negligible. Bound aflibercept might dominantly eliminate by internalization (k int >>k off ). This assumption is reasonable and consistent with the study of Eppler et al on the development of a TMDD model for recombinant human VEGF (rhVEGF) after i.v. administration in patients with coronary artery disease [START_REF] Eppler | A target-mediated model to describe the pharmacokinetics and hemodynamic effects of recombinant human vascular endothelial growth factor in humans[END_REF]. VEGF binding to its high affinity receptors was concluded as an essentially irreversible process in vivo. Aflibercept might have the same properties because this drug is produced from two high affinity receptors of aflibercept (VEGFR1 and 2). During the model development, the general TMDD model with k on and k off was also implemented, assuming that the total target concentration was constant. Free VEGF was then the difference between total VEGF and VEGF in binding to aflibercept. This model showed that k on and k off could not be estimated separately with very small value of k off and its poor precision.

Removing k off resulted in the model developed by Epper et al [START_REF] Eppler | A target-mediated model to describe the pharmacokinetics and hemodynamic effects of recombinant human vascular endothelial growth factor in humans[END_REF], but we obtained unsatisfying fits for bound aflibercept, either for central or peripheral binding to VEGF. A plateau phase after reaching the peak, instead of decreasing phase was predicted for bound aflibercept and the predicted peak was lower. 20 Alternatively, we also evaluated the linear central or peripheral binding models with a first-order association rate constant (k on ) instead of V max and K m . These models did not describe the observed concentrations of bound aflibercept well, specifically the decreasing phase which illustrate the elimination of bound alibercept after reaching the peak was poorly approached. They were also found less appropriate by statistical criteria.

The nonlinear peripheral binding model was found to better describe the kinetics of drugs than the central binding one, which is consistent with the previous studies on the distribution of VEGF in the body [START_REF] Kut | Where is VEGF in the body? A meta-analysis of VEGF distribution in cancer[END_REF]. Large quantities of VEGF were reported in extravascular space of tumour and skeletal muscle [START_REF] Kut | Where is VEGF in the body? A meta-analysis of VEGF distribution in cancer[END_REF], suggesting an important source of endogenous VEGF in peripheral compartment and supporting our choice of final model. The plasma concentrations of bound aflibercept in the data were assumed to be the same as those in peripheral compartment, with rapid transfer between extravascular and plasma space.

In the final model, the clearance of bound aflibercept was found to be 6.3 times lower than that of free aflibercept from central compartment (0.14 L/day and 0.88 L/day, respectively). Both free and bound aflibercept eliminate quite slowly. The typical central volume of distribution of free aflibercept was 4.94 L, indicating that free aflibercept has a low level of tissue diffusion and it circulates mostly in the extravascular spaces. The volume of distribution of bound aflibercept was close to that of free aflibercept, which was observed in intermediate models and was fixed to the mean value of V p in the final model due to the problem of identifiability. In the prior non compartmental analysis of free aflibercept, the average clearance was 0.97 L/day and the steady state volume of distribution was 5.98 L. These values were similar to those estimated by the modelling approach. Compared to a similar antiangiogenic agent, bevacizumab, the central volume of distribution is close but the clearance of total aflibercept is nearly 4 times faster than that of total bevacizumab [START_REF] Lu | Clinical pharmacokinetics of bevacizumab in patients with solid tumors[END_REF]. From MM parameters, the production rate of VEGF and the irreversible binding constant were calculated to be 0.99 mg/day and 2.91 µg/mL or 19909 pM, respectively. The irreversible binding constant is much larger than the in vitro binding affinity K d observed for VEGF-A. Such a difference is due to the negligible value of k off in vivo and it means that the target degradation rate (k deg ) is significantly faster than the dissociation rate (k off ).

The advantage of this study is the availability of drug concentration data under free and bound forms since they were separately assayed, which is not the case for many antibody products.

Generally, only total drug, representing the sum of free and bound was assayed [START_REF] Lu | Clinical pharmacokinetics of bevacizumab in patients with solid tumors[END_REF][START_REF] Bhargava | A phase I and pharmacokinetic study of squalamine, a novel antiangiogenic agent, in patients with advanced cancers[END_REF][START_REF] Hayashi | A mechanism-based binding model for the population pharmacokinetics and pharmacodynamics of omalizumab[END_REF]. The pharmacokinetic analysis of total form could not provide a good explanation of the mechanism of action of the drug, including the binding to the target; neither was it able to well characterise its pharmacokinetic profile. For example, a simple model of two-compartment infusion with firstorder elimination was published for bevacizumab and squalamine, two antiangiogenic agents dosed in total form [START_REF] Lu | Clinical pharmacokinetics of bevacizumab in patients with solid tumors[END_REF][START_REF] Bhargava | A phase I and pharmacokinetic study of squalamine, a novel antiangiogenic agent, in patients with advanced cancers[END_REF]. The availability of complex data, for aflibercept, helps us characterise both the linear and nonlinear elimination pathways, as well as the complex internalization, and estimate the mechanistic parameters of TMDD system; this could not be achieved using a classical MM model of free drug alone. The concentrations of bound aflibercept could be considered as pharmacodynamic data information of this drug. Our model development was in agreement with the guidelines of modelling for drugs with TMDD properties as mentioned by Yan et al [START_REF] Yan | Selection between Michaelis-Menten and target-mediated drug disposition pharmacokinetic models[END_REF]. To our knowledge, the model that we developed is the first mechanism-based population pharmacokinetic model for an anti-VEGF drug.

In conclusion, the present PK model for aflibercept well characterises the underlying mechanism of disposition of aflibercept, where a saturable, high-affinity binding of the aflibercept to its pharmacologic target (VEGF) is responsible for the observed nonlinear pharmacokinetic 22 behaviour of the free drug. Although further studies will be needed to assess the influence of covariates because of the limited size of the present sample, this model helps to better understand the properties of aflibercept and provides a useful support for further studies in patients during the clinical development, in particular the determination of therapeutic doses using bound aflibercept concentrations as a marker of VEGF inhibition.

Competing interests

Hoai-Thu Thai received a research grant from Sanofi-aventis to perform this research and is now doing a PhD funded on a research grant from Sanofi-aventis. Emmanuelle Comets acted as a consultant for Sanofi-aventis in 2009-2010. France Mentré is the director of Inserm UMR 738 which receives a research grant from Sanofi-aventis. Christine Veyrat-Follet and Gerard Sanderink are employees of Sanofi-aventis. Nicole Vivier and Catherine Dubruc are former employees of Sanofi-aventis.

APPENDIX: Model development

General peripheral model of target-mediated drug disposition [START_REF] Mager | General pharmacokinetic model for drugs exhibiting target-mediated drug disposition[END_REF][START_REF] Mager | Target-mediated drug disposition and dynamics[END_REF]:

k el k int Drug k tp k pt C p , V p C t , V t R, V R RC (C b ,V b ) + k syn k deg k on k off (A1)
where C p , C t , C b , R: the concentrations of free drug in central and peripheral compartment, complex (bound drug) and target; V p , V t ,: volume of distribution of free drug in central and peripheral compartment; V R , V b : volume of distribution of target and bound drug in peripheral compartment. k el : elimination rate of free drug from central compartment; k on , k off : association and dissociation rate constant; k int : internalization rate constant of complex; k syn , k deg : target production and degradation rate.

V V C k R C k R k k dt dR C k C k V V R C k dt dC V V C k V V R C k C k V V C k dt dC V V C k C k k Input dt dC . . . . . . . . . . . . . . . . . . . . ). ( deg int + - - = - - = + - - = + + - =
To write this general model, we did not assume that the drug, the target and the bound drug have the same volume of distribution, and we assume that the amount of bound drug created during a unit of time, produced in the compartment of distribution of the target, is k on .R.C t .V R [START_REF] Mager | General pharmacokinetic model for drugs exhibiting target-mediated drug disposition[END_REF][START_REF] Mager | Target-mediated drug disposition and dynamics[END_REF].

When the free drug concentrations are much higher than target concentrations (C>>R) and their binding results in a fully saturated target, the full TMDD can be replaced by the Michaelis-Menten-type (MM) approximations of TMDD model following the suggestions of Gibiansky et al [START_REF] Gibiansky | Approximations of the target-mediated drug disposition model and identifiability of model parameters[END_REF][START_REF] Gibiansky | Target-Mediated Drug Disposition: New Derivation of the Michaelis-Menten Model, and Why It Is Often Sufficient for Description of Drugs with TMDD[END_REF]. We explain these approximations for peripheral TMDD models.

Approximation 1:

The first MM approximation of TMDD model was proposed in case of reversible binding, basing on the assumptions that the drug-target complex is in a quasi-steady-state or the derivative of the complex concentration, dC b /dt, is zero [START_REF] Gibiansky | Approximations of the target-mediated drug disposition model and identifiability of model parameters[END_REF]. V R and V b were assumed to be equal to V t .

C K C R C K k k k C R C C k k R C k + = = + = = + - . . 0 ). ( . . int int t ss t tot t tp t p p pt t C K C R k C k V V C k dt dC + - - = . . . . . int V max = R tot . k int ; K m =K ss =(k int +k off )/k on (A6)
If, in addition, the total target concentration R tot can be considered constant, V max is a constant parameter of the system. The TMDD equation system then results in two following equations:

(A7) Approximation 2:
When the binding is irreversible, the dissociation binding rate constant k off =0 then the TMDD equations simplify:

(A8)
The second MM approximation of TMDD models was proposed in the case of irreversible binding. It is based on the assumption that the target is in a quasi-steady-state so that the derivative of target concentration, dR/dt, is zero [START_REF] Gibiansky | Target-Mediated Drug Disposition: New Derivation of the Michaelis-Menten Model, and Why It Is Often Sufficient for Description of Drugs with TMDD[END_REF].

(

C K C V C k V V C k dt dC V V C k C k k Input dt dC + - - = + + - = . . . . . . ). ( max t IB IB t on syn t on syn C K K R C k k k R R C k R k k + = + = = - - . . 0 . . . 0 deg deg R C k R k k dt dR C k V V R C k dt dC V V R C k C k V V C k dt dC V V C k C k k
- - = - = - - = + + - =
where K IB is the irreversible binding constant (K IB = k deg /k on ) and R 0 is the target concentration at baseline (R 0 =k syn /k deg ), which can be obtained assuming steady-state in the differential equation for R in the absence of drug.

Using this expression for R and noticing that R 0 .K IB .k on =k syn , we can simplify the equation for dC t /dt to: (A11) Aslo, generally, only free drug is measured, so that in their model, Gibiansky et al did not take into account the evolution of the complex (dC b /dt) [START_REF] Gibiansky | Target-Mediated Drug Disposition: New Derivation of the Michaelis-Menten Model, and Why It Is Often Sufficient for Description of Drugs with TMDD[END_REF]. The TMDD equation system therefore includes only two equations: (A12) Application to our study:

In our case, the free aflibercept concentrations were much higher than the target concentrations.

We also considered the evolution of bound aflibercept since the concentrations of bound aflibercept were available and could not be considered negligible in the system. In addition, the dissociation rate constant was found to be very small during model development so that irreversible binding and the negligible change in target concentrations (dR/dt) could be assumed, which is similar to the assumptions of the second MM approximation described above. Thus, our differential equation system is (A12) to which we add the equation for dC b /dt from (A8) where

k off =0. t IB t R syn t t tp t p p pt t C K C V k V C k V V C k dt dC + - - = . . . 1
C K C V k V C k V V C k dt dC V V C k C k k Input dt dC + - - = + + - = . . . 1
C k C K C V V dt dC C K C V V C k V V C k dt dC V V C k C k k Input dt dC . . 1 . 1 . . . . . ). ( int max max - + ⋅ = + ⋅ - - = + + - = b t IB t R syn b b t IB t R
C k C K C V k V dt dC C K C V k V C k V V C k dt dC V V C k C k k Input dt dC . . . 1 . . . 1 . . . . . ). ( int - + = + - - = + + - = t IB IB C K K R R + = .
Tables Table 1-

Observed concentrations of free aflibercept (µg/mL) and bound aflibercept (µg.eq/mL) versus time.

Figure 2 Proposed structural models for free and bound aflibercept (2 compartments for free aflibercept, 1 compartment for bound aflibercept) with binding to VEGF occurring in the peripheral compartment.

  The goodness-of-fit plots of the best joint model are shown in Figure3. The plots of observations versus population predictions and observations versus individual predictions indicated that the model adequately described the observations over the dose range. The normalised prediction distribution error plots, NPDE versus time and NPDE versus predictions, show a symmetric distribution around zero, except for the early times and the small concentrations of bound aflibercept. This bias was caused by omitting the NPDE corresponding to below the quantification limit (BQL) observations. For BQL data, the metric NPDE has not yet been developed. The total number of BQL concentrations of bound aflibercept predicted by the model was 251 (33.73%), 231 of which correspond to observed BQL data. This shows the good agreement between model prediction and observation for BQL data. Examples of individual fits taken from two studies with three doses are shown in Figure 4. The model describes the observed concentrations of free and bound aflibercept for all subjects quite well.

  Recently, a new derivation of MM approximation of TMDD, called the irreversible binding MM model has been proposed by Gibiansky et al when the dissociation rate constant is negligible and the free target concentration is much smaller than the free drug concentration [19]. The developed model of aflibercept is therefore close to this new approximation of TMDD model. The only difference is that an extra differential equation was added to describe the evolution of bound drug. The MM parameters in the final model represent the combinations of the TMDD model parameters: V max =k syn .V R and K m =K IB =k deg /k on , where V R is the volume of distribution of target and K IB , k syn , k deg are the irreversible binding constant, target production rate and target degradation rate, respectively. Although these estimates are not sufficient for the complete description of the full TMDD model, they can provide useful information about the underlying kinetics of drug and target.

  where R tot is the total concentration of target (R tot =R+C b )We can rearrange the equation for dC t /dt, by expressing C b as a function of C t and using the expression of K ss , yielding the following equaation:(A5)This new equation dC t /dt therefore includes a Michaelis-Menten-type elimination with

Figure 3 Diagnostic

 3 Figure 3

Figure 4 Examples

 4 Figure 4

Figure 5 VisualFiguresFigure 1 Time

 51 Figure 5

  . . .

	p						t	t
			el	pt	p	tp
							p
	t	p	p			syn	R	t
	pt		tp	t
		t			t	IB	t

  max =k syn .V R =A syn and K m = K IB For parameter estimation, the micro constants (k el , k pt , k tp ) can be replaced by the macro constants: Q=k tp .V t =k pt .V p and CL=k el .V p ; where CL is the clearance of free drug from central compartment and Q is the intercompartment clearance of free drug.

	Substituting the expression for						derived from the second assumption of TMDD
	model, we find the following expression:
													(A13)
	Again, we notice a Michaelis-Menten elimination for equation dC t /dt, which enters the equation of
	dC b /dt as a saturable input.
	This equation system can be written with MM parameters, which represents our final model:
	p			(	el		). pt	p	tp	.	t	p .	t
	t	. pt	. p	t	p	. tp	t	t		m		t	t	(A14)
	b				t
	where V b t m b

Parameter estimates for the best free and bound model

  Clearance of free aflibercept from central compartment (CL=k el .V p ) Q: Intercompartment clearance of free aflibercept (Q=k tp .V t =k pt .V p )

	Fixed effects	Interindividual variability	Residual variability
	Parameter	Estimate (RSE%)	ω (%) (RSE%)		σ a (µg/mL) (RSE%)	σ p (%) (RSE%)
	CL (L/day)	0.88 (4.0)	28.0 (10)	Free aflibercerpt	0.05 (9.0)	17.1(3.0)
	V p (L)	4.94 (4.0)	27.3 (10)	Bound aflibercept	-	12.6 (4.0)
	Q (L/day)	1.39 (9.0)	49.8 (14)
	V t (L)	2.33 (7.0)	39.8 (14)
	V b (L)	4.94 (=V p )	-
	V max (mg/day)	0.99 (5.0)	13.6 (17)
	K m (µg/mL) k int (day -1 )	2.91 (11) 0.028 (5.0)	45.6 (15 -

CL:

Acknowledgments

The authors thank Global Metabolism and Pharmacokinetics Department, Sanofi-aventis, Paris which supported Hoai-Thu Thai by a research grant during this work.

List of abbreviations


Christophe Collewet

Eric Marchand

Photometric visual servoing

This paper proposes a new way to achieve robotic tasks by 2D visual servoing. Indeed, instead of using classical geometric features such as points, straight lines, pose or a homography, as is usually done, the luminance of all pixels in the image is here considered. The main advantage of this new approach is that it does not require any tracking or matching process. The key point of our approach relies on the analytic computation of the interaction matrix. This computation is based either on a temporal luminance constancy hypothesis or on a reflection model so that complex illumination changes can be considered. Experimental results on positioning and tracking tasks validate the proposed approach and show its robustness to approximated depths, low textured objects, partial occlusions and specular scenes. They also showed that luminance leads to lower positioning errors than a classical visual servoing based on 2D geometric visual features.

I. INTRODUCTION

Visual servoing consists aims at controlling the motions of a robot by using data provided by a vision sensor [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]. More precisely, to achieve a visual servoing task, a set of visual features has to be selected from the image allowing to control the desired degrees of freedom (d.o.f). A control law is then designed so that these visual features s reach desired values s * . The control principle is thus to regulate the error vector e = s-s * to zero. To build the control law, the knowledge of the interaction matrix L s is usually required [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF].

Visual features are always designed from visual measurements m(p k ) (where p k is the camera pose at time k) which requires a robust extraction, matching (between m(p 0 ) and m(p * ), where p * is the desired camera pose) and realtime spatio-temporal tracking (between m(p k-1 ) and m(p k )). However, this process is a complex task, as testified by the abundant literature on the subject (see [START_REF] Marchand | Feature tracking for visual servoing purposes[END_REF] for a recent survey), and is considered as one of the bottleneck of the expansion of visual servoing. Thus several works focus to alleviate this problem. An interesting way to avoid any tracking process is to use non geometric visual measurements as in [START_REF] Malis | A unified approach to visual tracking and servoing[END_REF], [START_REF] Benhimane | Homography-based 2d visual tracking and servoing[END_REF] instead of geometric measurements as it is usually done. Of course, directly using non geometric visual features also avoids any tracking process. In that case, parameters of a 2D motion model have been used in [START_REF] Questa | Camera self orientation and docking maneuver using normal flow[END_REF]- [START_REF] Crétual | Visual servoing based on image motion[END_REF]. Nevertheless, such approaches require a complex image processing task.

In this paper we show that this tracking process can be totally removed and show that no other information than the image intensity (the pure luminance signal) needs to be considered to control the robot motion. Indeed, to achieve this goal we use as visual measurement and as visual feature the simplest that can be considered: the image intensity itself. We therefore call this new approach photometric visual servoing. In that case, the visual feature vector s is nothing but the image while s * is the desired image.

Considering the image intensity as a feature has been considered previously [START_REF] Nayar | Subspace methods for robot vision[END_REF], [START_REF] Deguchi | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF]. However, those works differ from our approach in two important points. First, they do not directly use the image intensity since an eigenspace decomposition is performed to reduce the dimensionality of image data. The control is then performed in the eigenspace and not directly based on the image intensity. Second, the interaction matrix related to the eigenspace is not computed analytically but learned during an off-line step. This learning process has two drawbacks: it has to be done for each new object and requires the acquisition of many images of the scene at various camera positions. Considering an analytical interaction matrix, as we propose avoids these issues. An interesting approach, which also directly considers the pixels intensity, has been recently proposed in [START_REF] Kallem | Kernelbased visual servoing[END_REF]. However, only the translations and the rotation around the optical axis have been considered (that is the 4 most simple d.o.f.) whereas, in our work, the 6 degrees of freedom are controlled. However, an image processing step is still required. Our approach does not require this step.

In this paper, we summarize several previous works. In [START_REF] Collewet | Visual servoing set free from image processing[END_REF], the analytic computation of the interaction matrix related to the luminance for a Lambertian scene is provided, only positioning tasks have been considered. In [START_REF] Collewet | Modeling complex luminance variations for target tracking[END_REF], this matrix has been computed considering a lighting source mounted on the camera and the use of the Blinn-Phong illumination model (a simplified model of the Phong model detailed in the next section), only tracking tasks have been considered. In [START_REF]Photometry-based visual servoing using light reflexion models[END_REF], the Phong model has been used, only positioning tasks have been considered. In addition, these works refer to [START_REF]Photometric visual servoing[END_REF] where details concerning analytic computations are given. Note that in [START_REF] Dame | Entropy-based visual servoing[END_REF], although this is also a direct visual servoing approach, the considered features used in the control law are very different. In this paper, we specifically focus on the way the visual servoing problem has been turned into an optimization problem. More precisely, we analytically analyse the cost function to minimize in order to derive an efficient control law. Moreover, additional experimental results than those described in our previous works are presented, as for example a comparison between classical 2D geometric visual features and the use of the luminance. We will show that by using luminance much lower positioning errors can be obtained.

The remainder of this paper os organized as follows, we first compute the interaction matrix of the luminance in Section II. Then, we reformulate the visual servoing problem into an optimization problem in Section III, and propose a new control law dedicated to our cost function to minimize. Section IV shows experimental results on various scenes for several tasks.

II. LUMINANCE AS A VISUAL FEATURE

The visual features considered in this paper are the luminance of each point of the image. More precisely, we consider as visual features the luminance I x at a constant pixel location x = (x, y) for all x belonging to the image domain and for a given pose p. Thus, we have

s(p) = I x (p) = (I 1• , I 2• , • • • , I N • ) (1)
where

I i• is the i-th line of the image. I x (p) is then a vector of size k = N × M where N × M is the size of the image.
As mentioned in the introduction, an estimation of the interaction matrix is required to control the robot motion. In our case, we are looking for the interaction matrix L Ix related to the luminance I x (t) at time t, that is

İx = L Ix v (2)
with v = (v, ω) where v is the linear camera velocity and ω its angular velocity.

Let us consider a particular point P (t) belonging to the scene which projects into the camera plane at the point p(t). P (t) is time varying either because the camera is moving with respect to the scene or because the scene is moving itself with respect to the camera. Let us note that, unless explicitly stated otherwise, all the quantities are expressed in the camera frame. The computation of the interaction matrix (2) requires to write the total derivative of the luminance I(p(t), t) in p at time t

İ(p(t), t) = ∇I(p(t), t) ⊤ ṗ(t) + ∂I(p(t), t) ∂t . (3)
However, considering that, at time t, the normalized coordinates of p(t) coincide with x, (3) becomes

İ(p(t), t) = ∇I x (t) ⊤ ẋ + İx (t) (4)
with ∇I x (t) the spatial gradient of I x (t) and ẋ the 2D velocity of p(t). Therefore, to explicitly compute the interaction matrix L Ix , an illumination model is required to estimate İ(p(t), t).

The simplest one is, of course, the one that lied on the temporal luminance constancy hypothesis [START_REF] Horn | Determining optical flow[END_REF], as it is the case in most of computer vision applications. In that case, we simply have İ(p(t), t) = 0 and it becomes straightforward to derive the interaction matrix from (4) and (2) (see [START_REF] Collewet | Visual servoing set free from image processing[END_REF] for further details). In that case, we obtain

L Ix = -∇I ⊤ L x (5)
where the interaction matrix L x related to x (i.e. such that ẋ = L x v) has been introduced

L x = -1/Z 0 x/Z xy -(1 + x 2 ) y 0 -1/Z y/Z 1 + y 2 -xy -x .
(6) Of course, because of the temporal luminance constancy hypothesis, ( 5) is only valid for Lambertian scenes, that is for surfaces reflecting the light with the same intensity in each direction. Besides, ( 5) is only valid when the lighting source is not moving with respect to the scene. In fact, it is well known that the temporal luminance constancy hypothesis can be easily violated [START_REF] Verri | Motion field and optical flow: qualitative properties[END_REF]. Therefore, to derive the interaction matrix, we have to consider a more realistic reflection model than the Lambert's one 1 . In this paper, we derive the interaction matrix from the Phong illumination model [START_REF] Phong | Illumination for computer generated pictures[END_REF]. This model is not based on physical laws, but comes from the computer graphics community. Although empirical, it is widely used thanks to its simplicity, and also because it is appropriate for various types of materials, whether they are rough or smooth.

According to the Phong model (see Fig. 1), the intensity I(p(t), t) at point p and at time t writes as follows

I(p(t), t) = K s cos η α + K d cos θ + K a . (7)
This relation is composed of a diffuse, a specular and an ambient component and assumes a point light source. The scalar K s describes the specular component of the lighting; K d describes the weight of the diffuse term which depends on the albedo in P (t); K a is the intensity of ambient lighting in P (t). Note that K s , K d and K a depend on P (t). θ is the angle between the normal to the surface n in P (t) and the direction of the light source l; α is the angle between r (which is l mirrored about n) and the viewing direction d; r can be seen as the direction due to a pure specular object. The parameter η allows to model the width of the specular lobe around r, this scalar varies as the inverse of the roughness of the material. Considering that r, d and l are normalized, we can rewrite [START_REF] Santos-Victor | Visual behaviors for docking[END_REF] as

I(p(t), t) = K s u 1 η + K d u 2 + K a (8)
where u 1 = r ⊤ d while we have u 2 = n ⊤ l. Note that these vectors are easy to compute, since we have d = -x x and r = 2u 2 nl with x = (x, y, 1).

Consequently, from [START_REF] Crétual | Visual servoing based on image motion[END_REF], it becomes easy to compute İ(p(t), t) involved in (4) (we assume here that the scene is only constituted by one material)

İ(p(t), t) = ηK s u η-1 1 u1 + K d u2 (9)
that leads to a general formulation of the optical flow constraint equation [START_REF] Horn | Determining optical flow[END_REF] considering the Phong illumination model

∇I ⊤ x L x v + İx = ηK s u η-1 1 u1 + K d u2 . (10)
Thereafter, by explicitly computing the total time derivative of u 1 and u 2 and writing u1 as u1 = L ⊤ 1 v and u2 as u2 = L ⊤ 2 v where L 1 and L 2 are 6-dimensional vectors, we obtain the interaction matrix related to the intensity at pixel x in the general case2

L Ix = -∇I ⊤ L x + ηK s u 1 η-1 L ⊤ 1 + K d L ⊤ 2 . (11)
To compute the vectors L 1 and L 2 involved in [START_REF] Kallem | Kernelbased visual servoing[END_REF] we have to explicitly express u1 and u2 . However, to do that, we have to assume some hypothesis about how n and l move with respect to the camera. Various cases have been studied in [START_REF]Photometric visual servoing[END_REF]. Due to lack of place, we report here only the case of an eye-in-hand robot system where the light source is mounted on the camera and where the interaction matrix is computed at the desired position. It is a very classical way to proceed [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]. Indeed, it avoids to compute on-line 3D information like the depths for example. We consider here this case. More precisely, we consider that, at the desired position the depth of all the points where the luminance is measured are equal to a constant value Z * . That means that we consider that the object is planar and that the camera and the object planes are parallel at this position. Moreover, we consider a directional lighting source. In these conditions, all computations done (see [START_REF]Photometric visual servoing[END_REF]), the interaction matrix related to the luminance writes simply as

L Ix = ηK s u 1 η-1 x x Z y Z - x 2 + y 2 Z y -x 0 -∇I ⊤ L x . (12)
where Z = Z * x 2 .

Note that this matrix requires the computation of ∇I, let us point out that it is the only image processing step necessary to implement our method.

On the other hand, even if the analytical computation of the vectors L 1 and L 2 is not straightforward in the general case, their final expression is very simple and easy to compute in this particular case.

III. VISUAL SERVOING CONTROL LAW

The interaction matrix associated to the luminance being known, the control law can now be derived. For that, we have turned the visual servoing problem into an optimization problem where the goal was to minimize the following cost function [START_REF] Collewet | Visual servoing set free from image processing[END_REF] C

(p) = 1 2 e 2 (13)
where e = I x (p) -I x (p * ). However, it is well known that in visual servoing, some image motions due to particular camera motions are not observable. Of course, to derive an efficient control law such camera motions must be avoided.

In the next section we will first prove that such motions exist, compute them, and then propose a control law that will ensure a high decrease of the cost function.

At the desired position, thanks to a first order Taylor series expansion of the visual features I x (p) around p * , an approximation of the cost function in a neighborhood of p * can be obtained (see [START_REF] Collewet | Visual servoing set free from image processing[END_REF], [START_REF]Photometric visual servoing[END_REF] for more details):

C(p) = 1 2 (v∆t) ⊤ H * (v∆t). (14)
where

H * = L I * x ⊤ L I * x . ( 15
)
is the Hessian matrix at p * . Since L I * x is analytically known, H * is also known and ( 14) can be easily evaluated. However, we will consider here only the case where the temporal luminance constancy hypothesis is valid and, first, when the camera and the object planes are parallel. In this case, we will denotes H * by H * . We will consider the more complex case when these planes are not parallel afterwards.

By considering the relation between the normalized coordinates x and their pixel value u = (u, v), a line of L I *

x writes simply at first order in h

L I * x = ∇I x /Z * , ∇I y /Z * , -h(m∇I x + n∇I y )/Z * , -∇I y , ∇I x , -(n∇I x + m∇I y )h (16)
where x = m h and y = n h have been substituted in ( 5) with m = u -u 0 , n = v -v 0 where (u 0 , v 0 ) is the principal point of the camera and h = 1/F with F = f /µ, f being the focal length and µ the size of a pixel (supposed to be square). Note that since F is a high value, the first order Taylor series expansion ( 16) is valid. From ( 16), H * is easily obtained

H * =                      h 11 Z * 2 h 12 Z * 2 h h 13 Z * 2 - h 12 Z * h 11 Z * h h 16 Z * h 12 Z * 2 h 22 Z * 2 h h 23 Z * 2 - h 22 Z * h 12 Z * h h 26 Z * h h 13 Z * 2 h h 23 Z * 2 0 -h h 23 Z * h h 13 Z * 0 - h 12 Z * - h 22 Z * -h h 23 Z * h 22 -h 12 -h h 26 h 11 Z * h 12 Z * h h 13 Z * -h 12 h 11 h h 16 h h 16 Z * h h 26 Z * 0 -h h 26 h h 16 0                      (17)
where the h ij are functions of the image gradients computed at each pixel that is not useful to detail.

From [START_REF] Horn | Determining optical flow[END_REF], it is easy to show that the rank of H * is 4 since the first line is obtained from the fifth line divided by Z * and the second line from the fourth divided by -Z * . Therefore, whatever the image content is, 0 is a double eigenvalue, and its associated eigenvectors denoted e 1 and e 2 are simply generated by the kernel of H * Ker H * = -Z * 0 0 0 1 0 , 0 Z * 0 1 0 0 .

That means that along any direction d = γ 1 e 1 + γ 2 e 2 (with γ 1 and γ 2 non null scalars), the approximated cost function [START_REF]Photometry-based visual servoing using light reflexion models[END_REF] does not vary and therefore that the true cost function [START_REF] Collewet | Modeling complex luminance variations for target tracking[END_REF] will slowly vary. Note that these motions coincide with what it is observed in practice, it is always possible to compensate a x (respectively y) axis translational motion with a y (respectively x) axis rotational motion to keep an image almost constant from a different point of view. In addition, note that d does not depend at all on the image content, moreover d is a constant value. The problem is now to find a direction that highly decreases the cost function. Since d is a constant value, we search a direction orthogonal to d . Such a direction can be simply given by ∇C (p) since near the desired pose p * we have

∇C (p) ⊤ d = (v∆t) ⊤ H * ⊤ d = 0 (18)
(we recall that d ∈ Ker H * ). Note that, in practice, this direction is also valid even quite far from p * as will be proved by the experimental results. Now we investigate the more complex case where the camera desired position is no more parallel to the object plane. In that case, the depths are given by 1/Z = ax + by + c and the matrix H * writes now as

H * = H * + H * ∠ where H * ∠ =                -2b 11 c h -2b 12 c h 0 b 12 h -b 11 h 0 -2b 12 c h -2b 22 c h 0 b 22 h -b 12 h 0 0 0 0 0 0 0 b 12 h b 22 h 0 0 0 0 -b 11 h -b 12 h 0 0 0 0 0 0 0 0 0 0                (19)
with b ij = aα ij + bβ ij where α ij and β ij are functions of all the image gradients computed at each pixels. As in the previous case where the camera desired position was parallel to the object plane, the rank of H * is still 4 and the same conclusions can be led: whatever the image content is, 0 is a double eigenvalue the associated eigenvectors of which are still generated by the kernel of

H * Ker H * = -1/c 0 µ 1 0 1 µ 2 , µ 3 1/c µ 4 1 µ 5 0
where the µ k (k ∈ {1, • • • , 5}) are functions of the h ij and the b ij . Those are not useful to detail.

Here again, that means that there are some directions d generated by Ker H * where the cost function (13) slowly varies. Moreover, as previously, ∇C (p) is an optimal direction since near p * , ∇C (p) and d are orthogonal

∇C (p) ⊤ d = (v∆t) ⊤ H * ⊤ d = 0. (20)

B. Design of the control law

We propose the following algorithm to reach the minimum of the cost function. The camera is first moved in the direction of ∇C to highly decrease the cost function and next, to a direction according to d to explore the remainder 2-dimensional subspace to reach its minimum. The first step can be easily done by using a steepest descent approach. However, if the direction of ∇C (p) is almost constant, its amplitude is not constant and may even vary very slowly in practice. To cope with this problem we propose to use the following control law

v = -v c ∇C(p init ) ∇C(p init ) . (21)
That is, a constant velocity with norm v c is applied in the steepest descent computed at the initial camera pose. Consequently, this first step behaves as an open-loop system. To turn into a closed-loop system, we first detect roughly the minimum along the direction of ∇C from a 3 rd order polynomial filtering of C(p) and then apply a control law formally equal to the one used in the Levenberg-Marquardt approach. We denote MLM this method in the remainder of the paper. The resulting control law is then given by (see [START_REF]Photometric visual servoing[END_REF])

v = -λ (H + µ diag(H)) -1 L Ix ⊤ (I x (p) -I x (p * )) (22)
where λ and µ are positive scalars, H is the Hessian matrix at p. More precisely, first, a high value for µ is used in ( 22) (typically µ = 1) to turn the control law into a steepest descent like approach 3 and to reach the minimum along the direction of ∇C. Once this minimum has been reached, a lower value is used (typically µ = 10 -2 ) to switch continuously to a Gauss-Newton (GN) control law (commonly used in visual servoing, see [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]) in order to explore the remainder 2-dimensional subspace generated by d and to reach the minimum of the cost function. This way to proceed ensure both a high convergence rate and a correct robot path.

Remarks about the stability: since redundant visual features (that is k > 6 a 6 d.o.f. robot) have been used, as it is also the case in classical visual servoing, only the local stability can be obtained (see [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF] for a proof). However, as pointed out in [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF], this domain can be quite large in practice. In addition, since we use redundant visual features, it is clear that the potential dimension of the null space can be high. However, that does not mean at all that all the motions that belongs to this null space are feasible, see for example [START_REF] Schramm | Extended 2d visual servoing[END_REF]. They use redundant visual features but prove that there are no local minima.

IV. EXPERIMENTAL RESULTS

In all the experiments reported here, the camera is mounted on a 6 degrees of freedom gantry robot. Control law is computed in real-time on a Core 2 Duo 3Gz PC running Linux. Images are acquired at 66Hz using an IEEE 1394 camera with a resolution of 320 × 2404 . The size of the vector I x is then 76800. Despite this size, the matrix L Ix can be computed at each iteration if needed.

A. Positioning tasks using the basic temporal luminance constancy model

We assume in this section that the temporal luminance constancy hypothesis is valid, i.e. we use the interaction matrix given in [START_REF] Questa | Camera self orientation and docking maneuver using normal flow[END_REF]. In order to make this assumption as valid as possible, a diffuse lighting has been used so that the luminance can be considered as constant wrt to the viewing direction. Moreover, the lighting is also motionless wrt the scene being observed. In this section, we will first compare the GN and MLM methods, then compare the photometric visual servoing with respect to a classical approach based on a matching and a tracking process and, finally, we will show that the photometric visual servoing is robust.

Comparison between the GN and the MLM method. The goal of the first experiment is to compare the control laws based on GN (the usual visual servoing control law) and MLM approaches when a planar object (a photo) is considered. The initial error pose was ∆p init = (5 cm, -23 cm, 5 cm, -12.5 • , -8.4 • , -15.5 • ) 5 . The desired pose was so that the object and CCD planes are parallel at Z = Z * = 80 cm. The interaction matrix has been computed at each iteration but assuming that all the depths are constant and equal to Z * (i.e. Eq. ( 5) with Z = Z * ), which is of course a coarse approximation.

Fig. 2a depicts the behavior of the cost functions using the GN method or the MLM method while Fig. 2b depicts the trajectories (expressed in the desired frame) when using either the GN or the MLM method. Fig. 2c and Fig. 2d depict respectively the camera velocity. The initial and final images are shown respectively, in Fig. 2e and Fig. 2f. First, as can be seen in Fig. 2a, both control laws converge since the cost functions vanish 6 . However, the time-to-convergence with the GN method is much higher than that of the MLM method. The trajectory when using the GN method is also shaky compared to the one of the MLM method (Fig. 2b). The velocity of the camera when using the MLM method is smoother than when using the GN method (Fig. 2d and Fig. 2c). This experiment clearly shows that the MLM method outperforms the GN one. Note that in both cases the positioning errors is very low, for the MLM method we obtained ∆p = (0.26 mm, 0.30 mm, 0.03 mm, 0.02 • , -0.02 • , 0.03 • ). It is very difficult to reach so low positioning errors when using geometric visual features as it is usually done. Indeed, this nice result is obtained because e is very sensitive to the pose p.

Comparison with a feature matching process. Considering such images, it is also possible to extract geometric features like SIFT or SURF and match them between current and desired images. In this experiment we use SURF features (with the OpenCV implementation) and use these points coordinates within a classical 2D visual servoing control law and within a robust one [START_REF] Comport | Statistically robust 2D visual servoing[END_REF]. The later allows to reject wrong matched features directly at the control level. 200 points have been extracted in the desired image and between 50 and 100 points are matched at each iteration and considered in the control law.

v x v y v z ω x ω y ω z (c) (d) (e) (f)
The goal is here to compare the precision of the positioning between these two classical approaches and the new proposed one. Table I shows the obtained results. As expected the positioning error is far lower using the new proposed approach. Translation error (that is t ) is only 0.07 mm (our robot precision is 0.1 mm) which has to be compared to the 0.87 mm and 1.28 mm using the two other methods. Indeed, in a classical approach, an extraction process obviously introduces errors in the features coordinates which implies imprecisions in the positioning task. Similar results have been obtained from other initial positions and other scenes.

Behavior with respect to partial occlusions. This experiment deals with partial occlusions. The desired object pose as well as the initial pose are unchanged. After having moved the camera to its initial position, an object has been added to the scene, so that the initial image is now the one shown in Fig. 3a and the desired image is still the one shown in Fig. 2f. Moreover, as seen in Fig. 3b and Fig. 3c, the object introduced in the scene is also moved by hand during the camera motion which highly increases the occluded surface. Despite that, the control law still converges (see Fig. 3f). Of course, since the desired image is not the true one, the error cannot vanish at the end of the motion (see Fig. 3f). Nevertheless, the final positioning error is not affected by the occlusions since we Robustness with respect to non planar scenes. The goal of this third experiment is to show the robustness of the control law wrt non planar scenes (see Fig. 4). This figure shows that large errors in the depth are introduced (the height of the castle tower is around 30 cm). The initial and desired poses are still unchanged. Fig. 5 depicts this experiment. Here again, the control law still converges (despite the fact that L Ix has been computed with a constant depth Z * = 80 cm) and the positioning error is still low since we have ∆p = (0.2 mm, -0.0 mm, 0.1 mm, -0.01 • , 0.00 • , 0.06 • ). Influence of the image content. The goal of these last set of experiments is to show that, even if the luminance is used as a visual feature, our approach does not depend too much on the texture of the scene being observed. Fig. 6 depicts the behavior of our algorithm for different planar objects (the initial as well as the desired pose is unchanged). As can be seen, the control law converges in each case, and even in the case of a low textured scene. Let us point out that similar positioning errors than for the first experiment have been obtained.

B. Tracking tasks

Our goal is now to perform a tracking task with respect to a moving object. That is, we have to maintain a rigid link between the target to track and the camera. Considering that the scene is moving, a specific illumination model has to be considered as explained in section II. A directional light-ring is located around the camera lens 7 . The scene is then no more illuminated by a diffuse lighting. The object is unchanged but it is attached to a motorized rail that allows to control its motion (see [START_REF]Photometric visual servoing[END_REF] for an object moved by hand). Although only one d.o.f of the object is controlled (with a motion that is completely unknown from the tracking process), the 6 d.o.f of the robot are controlled (the object velocity is 1 cm/s).

Since we have a constant target velocity, a simple integrator, as in [START_REF] Chaumette | Positioning of a robot with respect to an object, tracking it and estimating its velocity by visual servoing[END_REF], has been introduced in the control law to eliminate the steady state tracking error.

In this experiment, we use equation ( 12) to compute L Ix . In this relation, two parameters depending on the object surface are required: K s and η. However, in practice, a large domain of values for these parameters has led to good results. The same values (K s = 200 and η = 200) have been used for all our experiments and never been changed despite the fact that various material has been considered (glass, various plastics, metal, glossy and matt paper), see [START_REF] Collewet | Modeling complex luminance variations for target tracking[END_REF] and [START_REF]Photometry-based visual servoing using light reflexion models[END_REF] for other experiments using a complex illumination model. When the velocity is constant, the object is perfectly tracked, as can be seen on Figure 7a where e is depicted, despite the occurrence of a large specularity which shows the importance of using a complex illumination model. Let us note that without using this new model, experiments fail [START_REF]Photometric visual servoing[END_REF]. Error in the image remains small except when the object stops or accelerates (corresponding to the peaks in Figure 7a). For each pixel, except during accelerations and decelerations, |I x -I *

x | < 5. The camera velocity (see Fig. 7b) shows a pure motion along the x (± 1cm/s) axis that corresponds to the ground truth.

V. CONCLUSION

We have shown in this paper that it is possible to use directly the luminance of all the pixels in an image as visual features in visual servoing. To the best of our knowledge this is the first time that visual servoing has been handled without any complex image processing task (except the image spatial gradient required for the computation of the interaction matrix), nor learning step. Indeed, unlike classical visual servoing where geometrical features are used, using photometric visual servoing does not need any matching between the initial and desired features, nor between the current and the previous features. It is a very important issue when complex scenes have to be considered. Our approach has been validated on various scenes and various lightings for positioning or tracking tasks. Concerning positioning tasks, the positioning error is always very low and much lower than classical visual servoing based on 2D geometric visual features. Supplementary advantages are that our approach is not sensitive to partial occlusions and to coarse approximations of the depths required to compute the interaction matrix.

Fig. 1 .

 1 Fig.1. Quantities involved in the Phong illumination model[START_REF] Phong | Illumination for computer generated pictures[END_REF] (expressed here in the scene frame).

Fig. 2 .

 2 Fig. 2. First experiment. MLM vs. GN method (x axis in second). (a) Comparison of cost functions, (b) Comparison of camera trajectories, (c) Camera velocities (m/s or rad/s) for the GN method, (d) Camera velocities (m/s or rad/s) for the MLM method, (e) Initial image, (f) Final image.

Fig. 3 .

 3 Fig. 3. Second experiment. Partial occlusions (x axis in second). (a) Initial image, (b) Image at t ≈ 11 s, (c) Image at t ≈ 13 s (d) Final image, (e) Camera velocities (m/s or rad/s), (f) Cost function, (g) Ix -I * x at the initial position, (h) Ix -I * x at the end of the motion. have ∆p = (-0.1 mm, 2 mm, 0.3 mm, 0.13 • , 0.04 • , 0.07 • ). It is very similar with the previous experiments. This very nice behavior is due to the high redundancy of the visual features we use.

Fig. 4 .Fig. 5 .

 45 Fig. 4. The non planar scene used in Fig. 5 experiment.

Fig. 6 .

 6 Fig. 6. Fourth experiment. Same positioning task wrt to various objects. Objects considered (left column); Cost functions (right column) (x axis in second).

Fig. 7 .

 7 Fig. 7. Target tracking considering the complete interaction matrix that integrates specularity, diffuse and ambient terms (x axis in frame number). (a) Error Ix -I * x , (b) Camera velocity (m/s and radian/s), (c) Images at different time (left) and corresponding errors Ix -I * x (right).

TABLE I FINAL

 I POSITIONING ERROR: WE COMPARED THE PROPOSED APPROACH WITH APPROACHES BASED ON GEOMETRIC FEATURE EXTRACTION.

			t (in mm)		u θ (in degrees)
	Initial error	143.042 -177.517 12.496 -16.083 -10.139 -1.517
	Photometric VS	-0.027	0.042	-0.049	0.001	-0.001	-0.006
	SURF (2D VS)	0.486	0.558	-0.467	0.041	-0.046	-0.018
	SURF (2D robust VS)	1.072	-0.560	-0.418	-0.056	-0.027	-0.074

Indeed, the Lambert's model can only explain the behavior of non homogeneous opaque dielectric material[START_REF] Reichmann | Determination of absorption and scattering coefficients for non homogeneous media[END_REF]. It only describes a diffuse reflection component and does not take into account the viewing direction.

Note that we recover the interaction matrix -∇I ⊤ Lx associated to the intensity under temporal constancy (see (5)), i.e. in the Lambertian case (Ks = 0) and when u2 = 0 (the lighting direction is motionless with respect to the point P ).

More precisely, each component of the gradient is scaled according to the diagonal of the Hessian, which leads to larger displacements along the direction where the gradient is low.

Note that if a higher resolution is used, the computations time of the error vector and of the interaction matrix will be highly increased, it is thus better to decrease their size by decreasing the image resolution. There is no real advantage to use high resolution images for control issues.

The following notation has been used: ∆p = (t, uθ) where t describes the translation part of the homogeneous matrix related to the transformation from the current to the desired frame, while its rotation part is expressed under the form uθ where u represents the unit rotation axis vector and θ the rotation angle around this axis.

In fact, the cost functions do not exactly vanish, the mean error of the intensity levels is 2.2 (with a standard deviation of 0.4) at the end of the motion for the MLM method, this error is due to the acquisition noise and not due to the positioning error which, as we shall see, is very low.

The nature of the light is directional because LEDs with a small emission angle have been used.

ACKNOWLEDGEMENTS

The authors wish to thank Franc ¸ois Chaumette and Seth Hutchinson for their constructive comments.


M Kucmin

A Kajňaková

V A Feher

Starodub

A Radváková

O N Kazheva

A N Chekhlov

O A Dyachenko

M Kajňaková

A Feher

V A Starodub

M Kajėaková

the two-stack anionradical salt (N-Me-3,5-Di-Me-Py)(TCNQ) 2 (Py is pyridine)

Keywords: A. organic compounds, B. chemical synthesis, C. X-ray diffraction, D. magnetic properties

    Two-gap magnetic structure of the two-stack anion-radical salt (N-Me-3,5-Di-Me-Py)(TCNQ) (Py is pyridine)

Introduction

A c c e p t e d m a n u s c r i p t

Organic materials based on TCNQ are considered as a proper systems for the study of low-dimensional magnetic models due to their significant property -small modification of their crystal structure has a substantial impact on their physical properties. Recently we have studied magnetic and thermodynamic properties of samples (N-Me-2,5-di-Me-Pz)(TCNQ) 2 [START_REF] Radváková | Spin-Peierls Transition in Dimerized Stacks of Anion-Radical Salt (N-Me-2,5-di-Me-Pz)(TCNQ) 2 , (Pz is pyrazine)[END_REF] described with a dimerized Heisenberg spin 1/2 chain model and (N-Me-2,6-di-Me-Pz)(TCNQ) 2 [START_REF] Radváková | Quasi-two-dimensional ARS (N-Me-2,6-di-Me-Pz)(TCNQ) 2 , Pz is Pyrazine)[END_REF] described with a spin-ladder model. These compounds have very similar composition, in which only the position of the methyl group on the cation core is varying, but they differ significantly in the crystal structure, due to small changes in cation molecules.

It was theoretically predicted in the past that the low-dimensional systems may exhibit the spin-Peierls transition [3][START_REF] Pytte | Peierls instability in Heisenberg chains[END_REF][5]. This transition was experimentally observed for the first time in charge-transfer complexes (TTF)[M{S 2 C 2 (CF 3 ) 2 }] (TTF = tetrathiafulvalene) (M = Cu, Au) [6,7] and MEM-(TCNQ) 2 (MEM = methyl-ethylmorpholinium) [8]. Decades later, the spin-Peierls transition was observed also in inorganic compounds, firstly in CuGeO 3 [9], which is a main representant of spin-Peierls inorganic systems, and later also in (VO) 2 P 2 O 7 [START_REF] Garrett | Magnetic Excitations in the S = 1/2 Alternating Chain Compound (VO) 2 P 2 O 7[END_REF], which is a model system of interacting chains. Moreover, in the latter compound two energy gaps were observed [START_REF] Yamauchi | Two Gaps in (VO)2P2O7: Observation Using High-Field Magnetization and NMR[END_REF].

Two-gap systems are insufficiently described in the literature. The so-called double gap model was introduced for the first time in 1999 [START_REF] Garrett | Magnetic Excitations in the S = 1/2 Alternating Chain Compound (VO) 2 P 2 O 7[END_REF], for (VO) In the crystal structure of the studied compound (N-Me-3,5-di-Me-Py)(TCNQ) 2 two types of columns appear. This has led us to describe the magnetic properties with a model of dimerized chains with two gaps.

In this paper we present the synthesis and the crystal and molecular structure investigations of a (N-Me-3,5-di-Me-Py)(TCNQ) 2 single crystal. X-ray diffraction experiments were performed. Magnetic susceptibility was measured in the temperature range from 2 K to 300 K. Analysis of the low-temperature susceptibility revealed that two energy gaps occur in the magnetic excitations spectra and that the spin-Peierls transition could proceed in a rather unusual way.

Experimental

TCNQ was purified by the recrystallization from acetone and after that the gradient sublimation in argon atmosphere was used. Anion-radical salt (ARS) was gained by a standard method using reaction of TCNQ reduction by N-Me-3,5-di-Me-pyridine iodide, which we get from the reaction of methyliodide with 3,5-diMe-pyridine:

MeI + (3,5-diMePy) → [N-Me-3,5-diMePy] + I -; 1.5[N-Me-3,5-diMePy] + I -+ 2TCNQ → [N-Me-3,5-diMePy](TCNQ) 2 ↓ + 0.5[N-Me-3,5-diMePy](I 3 ).
One and half residue of iodide is necessary for binding of the forming iode. Precipitate of the salt was filtered and washed with diethyl ether and acetone mixture. In order to grow single crystals the recrystallization from acetonitrile/acetone mixture was performed. Black-violet needle crystals with length up to 2 mm were obtained.

A c c e p t e d m a n u s c r i p t

The X-ray study of the crystal structure of the ARS was performed at room temperature using an Enraf-Nonius CAD-4 (ω/2θ-scanning, Mo-Kα radiation, graphite monochromator) device. The crystal structure was solved by direct methods and subsequent Fourier synthesis using the SHELXS-97 software [START_REF] Sheldrick | SHELXS 97[END_REF]. Structure was refined by full-matrix least squares procedures using an anisotropic approximation for all non-hydrogen atoms within the SHELXL-97 [START_REF] Sheldrick | SHELXL 97[END_REF] software. Hydrogen atoms coordinates were calculated from the geometric condition. Nitrogen atom in each independent cation can be localized in one of three positions to which three methyl groups are bound. From these three possible positions we have calculated the one with maximal occupation. Main crystallographic data and some experimental characteristics are listed in table The temperature dependence of susceptibility was investigated in the temperature range from 2 K to 300 K in applied magnetic field 100 mT. Diamagnetic contribution to the susceptibility was estimated using Pascal's constants [START_REF] O'connor | [END_REF][START_REF]Software MaTra2.exe designed by Z.Trávníþek from Palacky University in Olomouc and J[END_REF] and then subtracted from the total susceptibility.

Results and Discussion

A c c e p t e d m a n u s c r i p t

Crystal structure of the salt is formed by four TCNQ anion-radicals A, B, C and D (Fig. 1a) and two crystallographicaly independent 3,5-dimethyl-N-methyl-pyridine cations E and F, situated in general positions (Fig. 1b).

Performed study showed that this salt has a layered structure (Fig. 2) in which layers formed by cations [3,5-di-Me-N-Me] + alternate along the a axis with layers consisting of TCNQ anion radicals. Anion layers separated by cation layers are inhomogeneous and consist, in turn, of the type A or B anion radicals (Fig. 3a), or the type C or D anion radicals (Fig. 3b).

Anions form layers along the c axis and irregular columns along the b axis (Fig. 3). The layers are formed by anions of only one type. In columns the anions alternate in

A c c e p t e d m a n u s c r i p t

Interplanar angles between planes laid through the carbon atoms of cations are 66.7° and 68.7° for the E and F type of cation, respectively.

In ARS with composition Ct + (TCNQ) 2 -(where Ct + is cation) tetramers usually form in the TCNQ columns due to the metal-insulator transition because when the conductivity band is filled up to 1/4 the period must, in conformity with Peierls theorem, increase fourfold [START_REF] Peierls | Quantum Theory of Solids[END_REF]. However, in our compound two types of structural dimers form -(AB) and (CD). For the interpretation of the physical properties of this ARS it is necessary to know the charge distribution on TCNQ molecules. This distribution was calculated by the method presented in [START_REF] Krycowski | Crystallographic studies and physicochemical properties of ʌ-electron compounds. III. Stabilization energy and the Kekule contributions derived from experimental bond length[END_REF]. From the data in Table 2a and 2b we calculated HOSE (Harmonic Oscillator Stabilization Energy) energies for quinoid and Kekule's TCNQ structures for A, B, C and D particles and their charges.

The partial charges are -0.10 e (A), -1.09 e (B), -0.52 e (C) and -0.89 e (D) (Table 2c).

With regard to the accuracy of this interpretation we suppose that structural dimers (AB) and (CD) carry the charge about -1 and consequently the spin 1/2 and therefore with decreasing temperature further ordering should occur in the (N-Me-3,5-Di-Me-Py)(TCNQ) 2 crystal, caused by the spin-Peierls transition, due to which spinless structure tetramers with the charge -2 form.

The temperature dependence of magnetic susceptibility of this ARS confirms stated presumption. The magnetic susceptibility Ȥ of the sample shows a broad maximum at T max = 26.5 K ± 1 K (Fig. 4). After substracting the diamagnetic contribution, approximated by Pascal's constants and found to be -13.587•10 -5 emu•mol -1 , the susceptibilty data were fitted by the Currie-Weiss law in the temperature region from 130 K to 300 K. We determined the value of the Curie constant C to be 0.299 emu•K•mol -1 and that of Weiss constant Ĭ to be -48.7 ± 2 K, confirming the

A c c e p t e d m a n u s c r i p t

presence of antiferromagnetic interactions. As can be seen in Fig. 4, this dependence is typical for quasi-one-dimensional systems with antiferromagnetic interaction.

The temperature dependence of the effective magnetic moment was fitted with formula (Fig. 5):

T P T P eff + ⋅ = 2 1 μ , ( 1
)
which yields the values P 1 = 1.6 ȝ B /K and P 2 = 28.3 K. The value of the effective moment 1.6 ȝ B corresponds to 0.89 unpaired electrons.

Considering two types of structure dimers in the structure of ARS we can expect the appearance of two gaps in the magnetic excitations spectra. Therefore the temperature dependence of magnetic susceptibility was described with the formula:

1 • © § Δ - + 1 • © § Δ - + = T T A T T A T C 2 2 1 1 exp exp * χ , (2)
where the first term describes the contribution from paramagnetic impurities (free fragments of magnetic structure) and the second and third terms describe the contributions from (AB) and (CD) type of dimers. This formula fits with parameters ɋ*

= 0.0186 emu•K•mol -1 , A 1 = 0.167 emu•K•mol -1 , A 2 = 0.127 emu•K•mol -1 , Δ 1 = 74.1 K
and Δ 2 = 26.6 K (Fig. 4) satisfactorily the temperature dependence of magnetic susceptibility. The amount of paramagnetic impurities was estimated from the value of the Curie constant and is equal to 0.9 %.

As we know, the value of the exchange integral decreases exponentially with the increase of the interplanar distance. Therefore such considerable difference between the gap widths for (AB) and (CD) type of dimers can be explained by different interplanar distances in dimers (3.19 Å in A-B dimers and 3.28 Å in C-D dimers).

A c c e p t e d m a n u s c r i p t 4. Conclusion

This study shows that in the presented ARS (N-Me-3,5-Di-Me-Py)(TCNQ) 2 an unusual dielectric state is observed in which Peierls dielectric phase forms in two grades. At room temperature two types of strucure dimers exist and doubling of the period is observed. At temperatures below 50 K further ordering occurs, with dimers coupling to tetramers. Two gaps in the magnetic excitations spectra correspond to two types of structure dimers. Nevertheless, to confirm our presumptions resulting from the crystal structure and magnetic susceptibility analyses, the measurements of magnetisation in high magnetic fields and EPR measurements should be performed.

A c c e p t e d m a n u s c r i p t

  in quite similar spin chains and the (VO) 2 P 2 O 7 can be considered as a typical gauge of unusual low-dimensional compounds.

  sequences (AB)(BA) or (CD)(DC) and structural dimers formed by anions of diverse types ((AB) and (CD)) can be distinguished, differing from each other only by the deviation angle. The dihedral angles, formed by planes of anions A-A, B-B, C-C and D-D, are 161.5°, 162.2°, 168.7° and 168.1°, respectively. Anions in dimers are practically parallel (dihedral angle between planes laid through all carbon atoms in TCNQ anion is equal to 0.4°). The distance between planes in (AB) dimers is 3.19 Å and in (CD) dimers is 3.28 Å. Anion-radicals TCNQ are characterized by considerable deviations of nitrogen atoms from the carbon plane, with values 0.07 ǖ, 0.10 ǖ, 0.16 ǖ and 0.15 ǖ for atoms N(1A), N(4B), N(4C) and N(4D), respectively. Cations are practically planar, with maximal deviation being only 0.01 ǖ and 0.02 ǖ (atoms N(1E) and ɋ(7F), respectively). Cation layers separated by anion layers are homogeneous and formed by cations of E or F type. Inside the layer the cations differ from each other in their orientation.

A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t Figure Captions Fig. 1 .

 Captions1 Figure Captions

Fig. 2 .Fig. 3 .

 23 Fig. 2. The crystal structure of (N-Me-3,5-Di-Me-Py)(TCNQ) 2 .

Fig. 4 .

 4 Fig. 4. Temperature dependence of magnetic susceptibility (open circles) fitted with

Fig. 5 .A c c e p t e d m a n u s c r i p t

 5 Fig. 5. Temperature dependence of effective magnetic moment (open circles) of ARS

A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t

  Figure 1a

Figure 2 A c c e p t e d m a n u s c r i p t Figure 3a A c c e p t e d m a n u s c r i p t Figure 3b A c c e p t e d m a n u s c r i p t

 23a3b Figure 2

Figure 4 A c c e p t e d m a n u s c r i p t

 4 Figure 4

Figure 5

 5 Figure 5

Table 1 .

 1 Basic crystalographic data Compound (N-Me-3,5-di-Me-Py)(TCNQ) 2

	Chemical formula	C 32 H 20 N 9
	Molecular mass	530.57
	Crystal system	Monoclinic
	Space group	Pɫ
	a/Å	27.578(4)
	b/Å	12.977(2)
	c/Å	7.796(1)
	α/°	90
	β/°	93.36(2)
	γ/°	90
	U/Å 3	2785.2(7)
	Z	4
	λ/Å	0.7107
	D calc , /g⋅ɫm -3	1.26
	μ/mm -1	0.080
	Number of measured reflections	5204
	Number of independent reflections	5204
	Number of observed reflections [F0 > 4σ(F0)] 2959
	Parameters refined	740
	(2θ) max , °	50.00
	Interval for h	-32 ≤ h ≤ 32

Table 2a .

 2a Bond distances in cations

		A	B	C	D
	N(1)-C(1)	1.134(6)	1.139(6)	1.151(7)	1.132(7)
	N(2)-C(2)	1.145(6)	1.140(6)	1.152(6)	1.136(6)
	N(3)-C(11)	1.137(6)	1.140(6)	1.155(7)	1.138(7)
	N(4)-C(12)	1.140(6)	1.141(6)	1.153(6)	1.128(7)
	C(1)-C(3))	1.439(6)	1.421(6)	1.411(7)	1.431(7)
	C(2)-C(3)	1.432(7)	1.414(6)	1.409(6)	1.429(7)
	C(3)-C(4)	1.377(8)	1.423(8)	1.400(8)	1.406(8)
	C(4)-C(9)	1.446(6)	1.409(6)	1.431(6)	1.419(6)
	C(4)-C(5)	1.448(6)	1.415(6)	1.435(6)	1.416(6)
	C(5)-C(6)	1.334(8)	1.381(8)	1.351(8)	1.370(8)
	C(6)-C(7)	1.449(6)	1.405(6)	1.437(6)	1.426(6)
	C(7)-C(10)	1.367(8)	1.420(8)	1.388(8)	1.392(8)
	C(7)-C(8)	1.445(6)	1.408(6)	1.425(6)	1.412(7)
	C(8)-C(9)	1.336(8)	1.365(8)	1.342(8)	1.366(8)
	C(10)-C(11)	1.436(6)	1.415(6)	1.405(7)	1.425(7)
	C(10)-C(12)	1.439(6)	1.416(6)	1.410(6)	1.432(7)
		A	B	C	D

Table 2b .

 2b Bond distances in anions

		E	F
	N(1)-C(1)	1.29(1)	1.33(1)
	N(2)-C(2)	1.37(1)	1.40(1)
	N(1)-C(6)	1.48(1)	1.49(1)
	C(1)-C(2)	1.35(1)	1.36(1)
	C(2)-C(3)	1.37(1)	1.35(1)
	C(2)-C(7)	1.45(1)	1.52(1)
	C(3)-C(4)	1.32(1)	1.41(1)
	C(4)-C(5)	1.35(1)	1.37(1)
	ý(4)-C(8)	1.52(1)	1.47(1)

Acknowledgements

This work was supported by the grants of Slovak Research and Development Agency under the contracts LPP-0102-06, No. APVV-VVCE-0058-07, No. APVV-0006-07 and VEGA -1/0159/09. The financial support of U.S.Steel -Košice is acknowledged.