Jhp Van Wel
email: jvanwel@maastrichtuniversity.nl

Kpc Kuypers

E L Theunissen

W M Bosker

K Bakker

J G Ramaekers

Van Wel

Blockade of 5-HT 2 receptor selectively prevents MDMA induced verbal memory impairment

Keywords: MDMA, ketanserin, pindolol, 5-HT 2, 5-HT 1, serotonin

MDMA or "ecstasy" has been associated with memory deficits during abstinence and intoxication. The human neuropharmacology of MDMA induced memory impairment is unknown. This study investigated the role of 5-HT 2A and 5-HT 1A receptors in MDMA induced memory impairment. Ketanserin is a 5-HT 2A receptor blocker and pindolol a 5-HT 1A receptor blocker. It was hypothesized that pretreatment with ketanserin and pindolol would protect against MDMA induced memory impairment. Subjects (N=17) participated in a double-blind, placebo controlled, withinsubject design involving 6 experimental conditions consisting of pretreatment (T1) and treatment (T2). T1 preceded T2 by 30 minutes. T1-T2 combinations were: placebo-placebo, pindolol 20mg-placebo, ketanserin 50mg-placebo, placebo-MDMA 75mg, pindolol 20mg-MDMA 75mg and ketanserin 50mg-MDMA 75 mg. Memory function was assessed at Tmax of MDMA by means of a word learning task, a spatial memory task and a prospective memory task.

MDMA significantly impaired performance in all memory tasks. Pretreatment with a 5-HT 2A receptor blocker selectively interacted with subsequent MDMA treatment and prevented MDMA induced impairment in the word learning task, but not in the spatial and prospective memory task. Pretreatment with a 5-HT 1A blocker did not affect MDMA induced memory impairment in any of the tasks.

Together, the results demonstrate that MDMA induced impairment of verbal memory as measured in the word learning task is mediated by 5-HT 2A receptor stimulation.

Introduction

3,4-Methylenedioxymethamphetamine (MDMA) is the primary psychoactive constituent in the popular party drug Ecstasy. Use of MDMA has consistently been associated with learning and verbal memory deficits in recreational and abstinent users (reviews: [START_REF] Cole | The pre-clinical behavioural pharmacology of 3,4methylenedioxymethamphetamine (MDMA)[END_REF][START_REF] Morgan | Ecstasy (MDMA): a review of its possible persistent psychological effects[END_REF][START_REF] Parrott | Human psychopharmacology of Ecstasy (MDMA): a review of 15 years of empirical research[END_REF]; the prototypical example being a reduction in performance on immediate and delayed word recall tasks. Although the existence of these persistent memory deficits in MDMA users has been well established, their neuropharmacology has been a matter of continuous debate. There is some evidence that these deficits are caused by MDMA induced neurotoxicity as indicated by depleted serotonin (5-HT) levels in MDMA users and by dose response relationships between the extent of exposure to MDMA and the magnitude of memory impairment. Yet, imaging studies demonstrated that memory impairment in ecstasy users was not associated with binding to cortical 5-HT transporters or duration of abstinence, which suggests that memory deficits in MDMA users occur independent of serotonergic neurotoxicity [START_REF] Reneman | Cortical serotonin transporter density and verbal memory in individuals who stopped using 3,4-methylenedioxymethamphetamine (MDMA or "ecstasy"): preliminary findings[END_REF].

A similar conclusion was drawn by Kuypers and[START_REF] Kuypers | Transient memory impairment after acute dose of 75mg 3.4-Methylene-dioxymethamphetamine[END_REF]2007 papers. They set out to inquire whether memory performance is directly related to 5-HT synaptic availability. They treated 18 recreational MDMA users with single doses of MDMA 75mg in a double blind, placebo controlled study and conducted both a verbal memory test and a spatial memory test during intoxication (i.e. 1-2 hrs post drug) and during withdrawal (i.e. 25-26 hrs post drug). Subjects displayed transient memory impairment during intoxication while they performed at placebo level in the withdrawal phase, when serotonin levels were depleted. These results suggest that memory performance during MDMA intoxication is not directly affected by changes in the availability of 5-HT in the synaptic cleft.

The apparent dissociation between general 5-HT levels and MDMA induced memory impairment raises the question whether MDMA exerts its detrimental effect on memory through other mechanisms than synaptic 5-HT-availability alone. Two interesting possibilities include direct and indirect activation of postsynaptic 5-HT 2A and 5-HT 1A receptors respectively. It is clear from animal studies that 5-HT 2A receptors mediate learning and memory processes [START_REF] Galizio | Effects of MDMA, methamphetamine and methylphenidate on repeated acquisition and performance in rats[END_REF][START_REF] Kay | Differential effects of MDMA and scopolamine on working versus reference memory in the radial arm maze task[END_REF][START_REF] Van Wel Trigo | MDMA modifies active avoidance learning and recall in mice[END_REF] although it is not completely obvious whether 5-HT 2A Van Wel receptor drugs achieve their facilitating and impairing effects through agonism, antagonism or inverse agonism [START_REF] Meneses | Involvement of 5-HT (2A/2B/2C) receptors on memory formation: simple agonism, antagonism, or inverse agonism?[END_REF]. Still, there appears to be some consensus that a 5-HT 2A blockade not only improves learning, but reverses poor memory consolidation conditions associated with, among others, dysfunctional serotonergic neurotransmission (Meneses, 2007a[START_REF] Meneses | Stimulation of 5-HT 1A , 5-HT 1B, 5-HT 2A/2C , 5-HT 3 and 5-HT 4 receptors or 5-HT uptake inhibition: short-and long-term memory[END_REF] as well. 5-HT 2A receptor agonists such as mCPP have been shown to decrease learning and memory in rats [START_REF] Meneses | Are 5-HT (1B/1D) and 5-HT (2A/2B/2C) receptors involved in learning and memory processes?[END_REF]. The same mechanism might be responsible for the decrement in learning and delayed recall following MDMA administration, as it has been previously demonstrated that MDMA possesses a moderate affinity for activating the 5-HT 2A receptor [START_REF] Sadzot | Hallucinogenic drug interactions at human brain 5-HT 2 receptors: implications for treating LSD-induced hallucinogenesis[END_REF].

A second candidate mechanism by which MDMA could cause its disadvantageous effects on memory is indirect activation of the post-synaptic 5-HT 1A receptor. The latter has also been shown to modulate memory performance in animals and humans [START_REF] Liechti | Psychological and physiological effects of MDMA ("ecstasy") after pretreatment with the 5-HT 2 antagonist Ketanserin in healthy humans[END_REF][START_REF] Meneses | Are 5-HT (1B/1D) and 5-HT (2A/2B/2C) receptors involved in learning and memory processes?[END_REF][START_REF] Yasuno | Inhibitory effect of hippocampal 5-HT 1A receptors on human explicit memory[END_REF]. PET studies have shown a negative correlation between memory function and 5-HT 1A receptor agonists binding in the hippocampus. More specifically, the 5-HT 1A agonist tandospirone dose dependently impaired performance in an immediate and delayed word recall task [START_REF] Yasuno | Inhibitory effect of hippocampal 5-HT 1A receptors on human explicit memory[END_REF]. The memory impairing effects of 5-HT 1A agonists thus appear identical to those produced by MDMA. It is not known whether MDMA also acts as a direct 5-HT 1A receptor agonist. However, it may achieve the same net result by indirectly stimulating postsynaptic 5-HT 1A receptors during the short-term increase in availability of synaptic 5-HT for the period of intoxication.

The present study was designed to investigate the role of 5-HT 2A and 5-HT 1A receptors in MDMA induced memory impairment. Recreational MDMA users performed a number of learning and memory tasks after receiving single doses of MDMA in a placebo-controlled study. In line with previous findings and assumptions, it is hypothesized that (1) an acute dose of MDMA will produce impairments on laboratory measures of learning and memory; (2) MDMA-induced memory impairment will be reversed by pre-treatment with ketanserin, a 5-HT 2A receptor antagonist, if impairment is related to direct or indirect stimulation of 5-HT 2A receptors; and (3) MDMA-induced memory impairment will be reversed by pretreatment with pindolol, a 5-HT 1A partial agonist, if impairment is related to direct or indirect stimulation of 5-HT 1A receptors.

Van Wel

Methods and Materials

Participants

A total of 17 healthy MDMA-users (9 male, 8 female), aged between 19 and 27 (mean (SD) 22.76 (2.75)) participated in the study. They were, with the exception of 2, all mild to moderate users of MDMA who reported to have taken the drug on 1 -65 separate occasions (mean 10.94) in the previous year. Overall, subjects reported to have taken MDMA 3 to 780 times in their lifetime (mean 72.4 times). When the heavy users were excluded from the analysis, the mean lifetime use of MDMA Subjects were paid for their participation in the study.

Treatments and procedures

Subjects participated in a double-blind, placebo controlled, within-subject design involving 6 experimental conditions consisting of pretreatment (T1) and treatment (T2). T1 preceded T2 by 30 minutes. T1-T2 combinations were: placebo-placebo, Van Wel pindolol 20mg-placebo, ketanserin 50mg-placebo, placebo-MDMA 75mg, pindolol 20mg-MDMA 75mg and ketanserin 50mg-MDMA 75 mg. Conditions were separated by minimum wash-out period of 7 days to avoid cross-condition contamination. The 75mg dose of MDMA was selected because it falls within the normal range of recreational use [START_REF] Gouzoulis-Mayfrank | Neurotoxicity of drugs of abuse--the case of methylenedioxyamphetamines (MDMA, ecstasy), and amphetamines[END_REF] and has been consistently

shown to impair memory performance and produce robust subjective mood changes in a number of previous studies from our lab [START_REF] Kuypers | Transient memory impairment after acute dose of 75mg 3.4-Methylene-dioxymethamphetamine[END_REF], 2006, 2007;[START_REF] Ramaekers | Involvement of inferior parietal lobules in prospective memory impairment during acute MDMA (ecstasy) intoxication: an event-related fMRI study[END_REF]. Doses of pindolol 20mg and ketanserin 50 mg represent regular therapeutic doses that block approximately 40% of 5HT 1A receptors and 91% of 5HT 2 receptors respectively [START_REF] Brogden | Ketanserin. A review of its pharmacodynamic and pharmacokinetic properties, and therapeutic potential in hypertension and peripheral vascular disease[END_REF][START_REF] Rabiner | Drug action at the 5-HT (1A) receptor in vivo: autoreceptor and postsynaptic receptor occupancy examined with PET and [carbonyl-(11)C]WAY-100635[END_REF][START_REF] Sharpley | Slow wave sleep in humans: role of 5-HT 2A and 5-HT 2C receptors[END_REF] Order of conditions was balanced over subjects and sessions. All subjects received a training session before onset of the experimental sessions in order to familiarize them with the tests and procedures. Treatments at T1 and T2 were administered using a double dummy technique to synchronize time of maximal drug concentrations (Tmax). Memory function was assessed by means of a number of memory tasks between 1.5-2 hrs after T2 (at Tmax), during which time subjects were allowed to watch TV or read a book. In addition, blood pressure and body temperature were assessed as safety measures. A summary of the test day is given in Table 2.

-insert table 2 -Subjects were asked to refrain from drugs at least a week before the start of the experiment and during the study. Subjects were not allowed to use alcohol on the day prior to an experimental session and were requested to arrive at experimental sessions well rested. Drug and alcohol screens were performed prior to experimental sessions upon arrival of the subject. Drug screens assessed for the presence of benzodiazepines, opiates, cocaine, marijuana, MDMA and (meth)amphetamine. Women were given a pregnancy test. Study treatments were only administered if subjects tested negative for drugs and pregnancy.

Van Wel

Memory tests

Three tests with demonstrated sensitivity [START_REF] Kuypers | Transient memory impairment after acute dose of 75mg 3.4-Methylene-dioxymethamphetamine[END_REF][START_REF] Ramaekers | Involvement of inferior parietal lobules in prospective memory impairment during acute MDMA (ecstasy) intoxication: an event-related fMRI study[END_REF] to the impairing effects of MDMA on memory were included in the memory task battery.

Word learning task

The Word Learning Task (WLT) is the Dutch language version of the standardized, clinically validated test for verbal memory [START_REF] Rey | The clinical examination in psychology / L'examen clinique en psychologie[END_REF]. The test begins with the sequential presentation of 30 monosyllabic common nouns. Each word is shown on the computer display for 2 seconds and the subject has to read them out loud. When the series ends, the subject is required to recall as many words as possible. The number correct is scored as the 1st trial score. Thereupon the same list is presented in the same manner on 2 successive occasions. Numbers correct are scored as before.

The main parameter is the sum score of correctly recalled words over 3 trials. After a 30 minute delay, the subject is shown a series of 30 words on the computer display, comprising of 15 words of the original set and 15 new words in random order. He/she was asked to respond as fast as possible to indicate whether the given word was one of the original set. The number of correct recognitions and the average speed (ms) of correct recognitions are recorded as the main dependent parameters [START_REF] Deelman | Measurements of aphasic disorders. A brief description of the SAN-battery[END_REF][START_REF] Schmitt | Non-serotonergic pharmacological profiles and associated cognitive effects of serotonin reuptake inhibitors[END_REF].

Spatial Memory Task

The spatial memory task is based on the spatial localisation task [START_REF] Vermeeren | Comparison of acute alprazolam (0.25, 0.50 and 1.0 mg) effects versus those of lorazepam 2 mg and placebo on memory in healthy volunteers using laboratory and telephone tests[END_REF]. It assesses short-term memory for non-verbal information. The subject is briefly shown a fixation point in the center of the computer display. Shortly thereafter, a target appears at a random location for 500ms. The subject has to memorize the position of the target and, using a computer mouse, relocate the cursor as accurately as possible over that position. The cursor appears either immediately upon target offset or after a delay of 2 or 4 sec. The subject depresses a button to indicate that the cursor is at the recalled position of the target. The test consists of 75 trials, divided equally among 3 response delays. The sequence of delays is random. Mean localisation error (mm) and average reaction time (ms) are the main parameters.

Van Wel

Prospective memory task

The prospective memory task is a new paradigm [START_REF] Ramaekers | Involvement of inferior parietal lobules in prospective memory impairment during acute MDMA (ecstasy) intoxication: an event-related fMRI study[END_REF]

Pharmacokinetics

Blood samples were collected before the start of the memory tasks at 1.5 hrs post T2.

Blood samples were centrifuged immediately and the serum was subsequently frozen at -20°C until analyses for pharmacokinetic assessments. MDMA, pindolol and ketanserin concentrations were determined using solid phase extraction and gas chromatography with mass spectrometric detection.

Van Wel

Statistics

The hypotheses that pretreatment with ketanserin or pindolol would interact with MDMA induced memory impairment was tested in 2 separate General Linear Model

Results

Seventeen complete data sets entered statistical analysis. Statistical analysis of memory data is described in detail below. A summary of main effects on measures of memory performance is given in Table 3.

Word learning task

MDMA significantly decreased immediate recall in both the ketanserin and pindolol comparisons (F

Prospective Memory Task

Van Wel MDMA significantly increased the number of prospective memory failures in the NoGo trials in one of the GLM comparisons (F 1,16 = 7.8; p=.013) and showed a trend in the other (F 1,16 = 4; p=.06). The factors Ketanserin and Pindolol or their interaction with MDMA did not affect prospective memory. Mean (se) prospective memory failures as a percentage of the total prospective memory trials in each treatment condition are given in Figure 1.

-----Insert Figure 1 3-------

Safety data

Safety data revealed a significant main effect of MDMA (F 1,16 = 24.5; p<.001) , pindolol (F 1,16 = 26.7; p<.001 and ketanserin (F 1,16 = 4.48; p<.05) on blood pressure.

Overall, mean blood pressure ranged from 73,6 to 87,3 mmHg over all treatments.

MDMA increased mean blood pressure by 7 mmH relative to placebo, whereas ketanserin and pindolol mildly reduced mean blood pressure relative to placebo.

However, ketanserin and pindolol did not interact with the effect of MDMA on blood pressure. Mean body temperature ranged from 36,2 -36,8 °C over all treatment conditions. Overall, MDMA increased body temperature (F 1,16 = 4.4; p<.05). Relative to placebo, body temperature increased by 0.2 °C after MDMA alone. Pretreatments did not affect body temperature or the effect of MDMA on body temperature.

Van Wel

Pharmacokinetics

Mean drug concentrations in all treatment conditions are shown in Table 4.

---Insert Table 4----

Discussion

The goal of the present study was to investigate the role of 5-HT 2A and 5-HT 1A receptors in MDMA induced memory impairment. Subjects were given single doses of MDMA combined with placebo, ketanserin or pindolol, after which several learning and memory tasks were performed.

Single doses of MDMA produced memory impairments in all memory tasks.

GLM analyses revealed that overall, MDMA decreased immediate recall in the WLT, increased the number of prospective memory failures in the Prospective Memory Task and increased localisation error in the Spatial Memory Task. These results are fully in line with previous studies that have reported identical memory impairments during intoxication with MDMA (Kuypers andRamaekers, 2005, 2007;[START_REF] Ramaekers | Involvement of inferior parietal lobules in prospective memory impairment during acute MDMA (ecstasy) intoxication: an event-related fMRI study[END_REF] as assessed in identical memory paradigms. These results basically indicate that the current memory paradigms were sensitive to the impairing effects of MDMA and should serve as reliable models to assess changes in MDMA induced memory impairments after pretreatments with 5-HT 1A and 5-HT 2A receptor blockers, the main goal of the present study.

The hypotheses that pretreatment with ketanserin or pindolol would interact with MDMA induced memory impairment was tested in 2 separate GLM analyses that assessed main effects of MDMA, Ketanserin and MDMA x Ketanserin and main effects of MDMA, Pindolol and MDMA x Pindolol respectively. The GLM models were chosen because they offer a robust and direct assessment of the interaction between pretreatments and MDMA in only 2 major analyses. This procedure was much preferred over the alternative approach of multiple comparisons between drugdrug and drug-placebo that can only provide indirect evidence of pretreatmenttreatment interactions. Consequently, the present study only conducted drug-placebo contrasts as supportive, secondary tests following a significant interaction between pretreatment and MDMA as evinced by GLM ANOVA.

Van Wel

GLM analyses revealed only one significant interaction between pretreatmenttreatment in the 3 memory tasks. In the WLT, pretreatment with ketanserin prevented impairment of word learning during treatment with MDMA. Separate drug-placebo contrasts also provided indirect and supportive evidence of the interaction between ketanserin and MDMA. T1-T2 combinations of placebo-MDMA and pindolol-MDMA both produced significant reductions in immediate recall. Drug-placebo contrasts also showed that performance in the word learning task after combination of ketanserin and MDMA did not differ from that during the placebo-placebo combination. Together these data strongly suggest that MDMA-induced impairments on verbal memory are in large part caused by direct or indirect stimulation of the 5-HT 2A receptor. The present finding appears in line with a previous report on 5-HT 2A receptors densities in recent MDMA and ex-MDMA users [START_REF] Reneman | The acute and chronic effects of MDMA ("ecstasy") on cortical 5-HT 2A receptors in rat and human brain[END_REF].

That study showed that 5-HT 2A receptor densities were significantly lower in recent Pretreatment with pindolol did not affect memory impairments in any of the memory tasks after treatment with MDMA. The interaction between pindolol and MDMA never reached statistical significance in any of the GLM analyses. These data suggest that 5-HT 1A receptors do not play a crucial role in the neuropharmacology of MDMA induced memory impairment. These finding confirm recent results from a study assessing the role of 5-HT 1A receptors on MDMA effects on visual-spatial working memory as measured with a CANTAB test battery [START_REF] Hasler | Investigation of serotonin-1A receptor function in the human psychopharmacology of MDMA[END_REF] in a placebo-controlled study. A single dose of MDMA significant impaired visual-spatial memory as compared to placebo. However, pretreatment with pindolol did not significantly alter MDMA induced impairment of memory. It should be noted however that Pindolol is not a very selective drug. It is basically a beta-blocker that Van Wel also blocks approximately 40% of 5HT 1A receptors in the brain [START_REF] Rabiner | Drug action at the 5-HT (1A) receptor in vivo: autoreceptor and postsynaptic receptor occupancy examined with PET and [carbonyl-(11)C]WAY-100635[END_REF].

Unfortunately, 5HT 1A ligands that selectively and fully block 5HT 1A receptors are presently not available. Nevertheless, recent molecular imaging studies also suggested that 5-HT 1A receptors may not be involved in human cognition at all [START_REF] Borg | Molecular imaging of the 5-HT 1A receptor in relation to human cognition[END_REF][START_REF] Borg | Search for correlations between serotonin 5-HT 1A receptor expression and cognitive functions--a strategy in translational psychopharmacology[END_REF].These studies failed to find any correlation between regional

[11C]WAY 100635 binding to 5-HT 1A receptors and cognitive performance in humans as measured with standard neuropsychological tests batteries [START_REF] Borg | Molecular imaging of the 5-HT 1A receptor in relation to human cognition[END_REF][START_REF] Borg | Search for correlations between serotonin 5-HT 1A receptor expression and cognitive functions--a strategy in translational psychopharmacology[END_REF]. Together, these results question the role of 5-HT 1A receptors in human memory, as it cannot explain the general memory deficit in MDMA users.

The present study did not provide clues to neuropharmacological mechanisms of MDMA induced impairments in every memory domain. The preventive effect of ketanserin on MDMA induced memory impairment was selective for verbal memory performance in the Word Learning Task. Ketanserin did not affect spatial memory and prospective memory impairments following MDMA. The selectivity of the interaction between ketanserin and MDMA on verbal memory was also supported by the lack of any interactions between pretreatment and MDMA on safety measures such as blood pressure and body temperature. The latter finding is in line with a previous study that demonstrated that co-administration of ketanserin did not alter MDMA effects on physiological measures (i.e. blood pressure, heart rate and body temperature) and subjective measures (e.g. well being and positive mood) [START_REF] Herin | Role of the serotonin 5-HT 2A receptor in the hyperlocomotive and hyperthermic effects of (+)-3,4-methylenedioxymethamphetamine[END_REF][START_REF] Liechti | Psychological and physiological effects of MDMA ("ecstasy") after pretreatment with the 5-HT 2 antagonist Ketanserin in healthy humans[END_REF].

The selective interaction between ketanserin and MDMA effects on verbal memory may indicate differences in neural memory networks and differences in distributions of 5-HT 1A and 5-HT 2A receptors within these networks. It is clear from the functional neuroimaging literature that frontal, parietal and occipitotemporal cortex are generally involved in a large variety of memory tasks, but that the pattern of brain recruitments within these networks may differ as a function of memory task or memory process [START_REF] Linden | The working memory networks of the human brain[END_REF]. Autoradiographic mapping of 5-HT 1 and 5-HT 2 receptors in rhesus monkeys demonstrated complementary patterns of distribution of 5-HT 1 and 5-HT 2 receptors in cortical areas of frontal, parietal and occipital lobes. 5-HT1 receptors were concentrated in layers I-III, whereas 5-HT 2 receptors were primarily concentrated in throughout layers III and IV [START_REF] Lidow | Quantitative autoradiographic mapping of serotonin 5-HT 1 and 5-HT 2 receptors and uptake sites in the neocortex of the rhesus monkey[END_REF].

Involvement of 5-HT 2A receptors in MDMA induced memory impairment thus may depend on task related network activations that occur within these layers.

Van Wel

Recent studies have also shown that a functional genetic polymorphism of the 5-HT 2A receptor (His452Tyr) has been associated with variations in memory performance, more specifically in hippocampus-dependent memory tasks [START_REF] Schott | Genetic variation of the serotonin 2A receptor affects hippocampal novelty processing in humans[END_REF][START_REF] Sigmund | Fine-mapping at the HTR 2A locus reveals multiple episodic memory-related variants[END_REF]. Future studies may focus on whether genetic variation of the 5-HT 2A receptors affects the magnitude of MDMA induced verbal memory impairment or determines the sensitivity of MDMA users for MDMA induced impairment of verbal memory.

In sum, the present study demonstrated that single doses of MDMA impaired memory in a range of tasks, and that pretreatment with ketanserin offered selective protection against MDMA induced verbal memory impairment. This demonstrates that MDMA induced verbal memory impairment is mediated by 5-HT 2A receptor stimulation.

Van Wel

  that was developed to examine prospective memory performance in an event based memory task. The foreground task consists of 100 successive presentations of a letter (A or B) in the center of a computer screen. Subjects are required to respond to each letter as quickly as possible by pressing one of two response buttons. One button is pressed to indicate that the letter "A" appeared and the other to indicate the letter "B". Both letters are presented equally often. Subjects are informed about the trial number by means of a trial counter that is always present in the left top corner of the screen. In addition, subjects are presented at irregular times with a future trial number in the right top corner of the display. Subjects are instructed to remember this future trial number and withhold from responding to the foreground task during the actual occurrence of the future trial. The memory set of trial numbers is dynamic and contains up to 3 future trial numbers. A trial number in the memory set is replaced by a novel future trial number whenever the actual trial number matches a future trail number in the set. Trials during which subjects are expected to respond are classified as Go trials. Trials during which subjects are instructed to withhold a response are classified as No Go trials (prospective memory trials). Time between presentation of a future trial number and the actual occurrence of the trial (i.e. memory delay) varies between 1, 2 and 3 minutes, equally divided over all No go trials. In total, the prospective memory task consists of 68 Go trials and 24 No go trials. Number of correct prospective memory recalls (i.e. number of correct response inhibitions) in No go trials is the primary dependent performance parameter.

  Figure 2. ----Insert Fig 2 andTable 3-------

  MDMA users and significantly higher in ex-MDMA users. The authors speculated that the latter was evidence of compensatory upregulation of post-synaptic 5-HT 2A receptors due to low synaptic 5-HT levels. The present study suggests that down regulation of 5-HT 2A receptors in recent MDMA users may result from high synaptic levels of 5-HT and 5-HT 2A receptor stimulation during MDMA intoxication. Studies showing involvement of 5-HT 2 receptors in acute MDMA intoxication as well as in persistent memory impairments in chronic MDMA users indicate that these phenomena may share related underlying mechanisms. Any definitive conclusion in that direction however remains premature since we do not know at present whether acute and chronic memory effects of MDMA are truly interconnected.

Figure 1 .

 1 Figure 1.Mean (SE) immediate recall and prospective memory recall in the word learning task and the prospective memory task for every treatment condition. Each treatment condition consisted of a pre-treatment (T1) and a treatment (T2). (PLA=placebo, PIN=pindolol, KET=ketanserin, MDMA=3,4-methylenedioxymethamphetamine, * = p<.05 as indicated by drug-placebo contrasts following a significant main interaction effect of MDMA x Ketanserin in the Word Learning Task).

Table / Figure legends Table 1 .

 /legends1 Subject demographics -drug history (lifetime use in number of times used).

		Mean	SD	Min	Max
	Age (years)	22,76	2,75	18,6	27,2
	Height (cm)	179,6	12,3	160,0	200,0
	Weight (kg)	72,9	14,7	56,0	110,0
	Alcohol	1168,5 794,1	260,0	2600,0
	Cannabis	554,4	1104,8 0,0	3650,0
	Ecstasy	72,4	194,7	3,0	780,0
	Amphetamines 4,1	12,3	0,0	49,0
	Cocaine	4,7	7,3	0,0	25,0
	LSD	4,4	12,0	0,0	50,0
	Mushrooms	0,2	0,5	0,0	2,0
	Other	0,5	1,7	0,0	7,0
	Van Wel

Table 4 .

 4 Dose and mean (SD) drug concentrations in each of the 6 treatment conditions. During pre-treatment subjects received either placebo, ketanserin 50mg or pindolol 20mg. Treatment consisted of either MDMA 75mg or placebo. Each treatment condition consisted of a pre-treatment (T1) and a treatment (T2).

	(PLA=placebo,	PIN=pindolol,	KET=ketanserin,	MDMA=3,4-
	methylenedioxymethamphetamine).
	T1	PLA	PIN	KET	PLA	PIN	KET
	T2	PLA	PLA	PLA	MDMA	MDMA	MDMA
	Dose (mg)	-	20	50	-	20	50
		-	-	-	75	75	75
	Cmean (ng/mL)	-	133 (80) 86 (42)	-	130 (53) 104 (41)
		-	-	-	157 (48) 156 (56) 164 (62)

Acknowledgements

The authors thank Dr S Toennes and Dr G Kauert (Institute of Forensic Toxicology, Goethe University Frankfurt) for establishing drug-serum concentrations and MSc Uschi Poesmans for her assistance in data collection.

Van Wel

Disclosure/Conflicts of interest

This research was funded by the Netherlands Organisation for Scientific Research (NWO), grant number: 400-05-096, awarded to JR. The authors declare that over the past 3 years they have also received funding for research not related to this paper. JR, ET and KK have received grants from pharmaceutical industries, the European Committee, and the Netherlands Organisation for Scientific Research (NWO).


M Sega
email: sega@science.unitn.it

E Autieri

F Pederiva

enhanced sampling of six-membered ring conformations

Keywords: Puckering, Sugars, Ring conformers, Free energy Puckering, Sugars, Ring conformers, Free energy

Pickett angles and Cremer-Pople coordinates as collective variables for the enhanced sampling of six-membered ring conformations.

Introduction

Cyclic structures are ubiquitous in nature, and appear mostly in non-planar, puckered conformations. The problem of an intuitive description of puckered ring conformers dates back to the work of Sachse [1], but it is only in 1975 that an exhaustive, mathematically correct description of a general N -membered ring puckering has been proposed [2]. The description of Cremer and Pople (CP) has from time to time been questioned because of some conflicts with the standard stereochemical interpretation of puckering [3], of its not straightforward way to determine atomic positions from puckering parameters [4][5][6], and for some difficulties in comparing with NMR experimental results [5,7]. Nevertheless, the CP approach has become the most widely employed way to describe puckered conformers. The calculation of puckering free energy landscapes, and the determination of conformer populations has been the subject of a series of recent computer simulation investigations [8][9][10][11][12], apparently stimulated by the possibilities offered by recent developments in accelerated dynamics methods. Among those works, Strauss-Pickett (SP) dihedrals [13][14][15] have been also employed, as an alternative to CP variables, to perform the exploration of the puckering free energy space.

In this work we address the advantages and disadvantages related to the choice of CP and SP coordinates, from the perspective of who is interested in computing the population of puckered conformers by means of computer simulations, in particular with the help of accelerated dynamics methods. (a) The Hammer-Aitoff equal-area projection [START_REF] Bugayevskiy | Map projections: a reference manual[END_REF] of the CP conformational globe. Conventional names and ring structures of some hexopyranose conformers are shown.

(b) Definition [15] of the out-ofplane dihedral angles (SP puckering representation).

Figure 1. Location of ring conformers on the CP globe and definition of SP dihedrals.

Puckering Coordinates and accelerated sampling methods

Puckering coordinates

The description of ring conformation can be either qualitative (that is, by assigning conventional names to special ring conformers) or quantitative (employing a parametric measure of the deviation from a reference shape). While IUPAC indications [16,17] realize the first of these two descriptions, no agreement on a standard quantitative measure has been reached so far. The first attempt to provide a set of parameters to describe ring conformers has been the observation of Kilpatric and coworkers [18,19], who employed for the first time the concept of a puckering amplitude and a pseudorotation angle to describe the puckered conformers of five-membered rings. A more systematic study of this problem came with the group theory analysis of Pickett and Strauss [20,21] and eventually with the first general definition of ring puckering coordinates of Cremer and Pople [2,6]. Both SP and CP approaches describe the conformation of a N -membered ring with N -3 parameters, expressed in terms of the elevations of the ring atoms from a reference plane (the mean ring plane for CP coordinates, and the planar ring conformation for SP coordinates). These two approaches, proven to be equivalent [22], represented the first mathematically consistent and quantitative measure of ring conformers for the puckering problem.

Cremer-Pople coordinates [2] are expressed in terms of the atomic elevations z j from the mean ring plane defined using the Cartesian coordinates of all ring atoms. In the case of six-membered rings the elevations from the mean plane,

z j = 1 √ 3 q 2 cos φ 2 + 2π(j -1) 3 + 1 √ 6 q 3 (-1) j-1 , (1)
depend on three puckering parameters, q 2 , φ 2 and q 3 . These coordinates are usually re-written as a set of spherical coordinates (Q, θ, φ), so that

   q 2 cos φ 2 = Q sin θ cos φ q 2 sin φ 2 = Q sin θ sin φ q 3 = Q cos θ, (2)
In this case the whole conformational space (see Fig. 1a) is spanned by amplitude, is related to the elevations by the simple relation

Q > 0, θ ∈ [0,
Q 2 ≡ N j=1 z 2 j = m q 2
m . Chair conformers are located on the north and south pole of the puckering sphere, while other conformers like boats, and skew boats are found along the pseudorotational path at the equatorial line.

The analysis of Pickett and Strauss, on the contrary, relied on the definition of the elevation Z j from a planar ring structure, and used group theory to identify the three parameters that can describe the puckering of a six-membered ring. Three out-of-plane dihedrals, α 1 , α 2 , α 3 (see Fig. 1) can be defined to be the three puckering parameters. For sugar rings, these dihedral angles are defined as:

α 1 =C 4 - O 5 -C 2 -C 1 -180 • , α 2 =O 5 -C 2 -C 4 -C 3 -180 • , α 3 =C 2 -C 4 -O 5 -C 5 -180 • ,
where the standard notation for ring atoms (C1 being the anomeric carbon) has been used. The SP description differs form the CP one only in the usage of the Z = 0 plane rather than the z = 0 mean plane, while the other the principles are the same [22]. These angles measure the orientation of tree "flaps" involving different atoms that can be identified on a six-membered rings (see Fig. 1b). Ring conformers are usually located in the [-90

• , 90 • ]×[-90 • , 90 • ]×[-90 • , 90 • ]
sub-volume. Common conformers for carbohydrates are listed in Tab.1. The three out-of-plane angles α i can also be rewritten [20] as

α i = R[cos Θ + 2 sin Θ(4πi/3 -2π/3 + Φ)], (3)
where again a spherical set of coordinates (R, Θ, Φ) is employed to describe the conformational phase space. It is worth remarking that this spherical coordinate set is different from the (Q, θ, φ) one, again because the two approaches differs in the definition of the reference plane.

Accelerated sampling methods

A large number of accelerated dynamics methods, such as metadynamics [START_REF] Laio | [END_REF], adaptive umbrella sampling [25], or local elevation [26] to cite but a few, have been proposed in the last decade to address the problem of investigating rare events (discussions regarding different enhanced sampling techniques can be found in Refs. [10,27]). The introduction of a time dependent biasing force which enforces the exploration of the complete phase space, at the price of loosing the physical meaning of microscopic trajectories, is the common trait of accelerated dynamics methods. The set of collective variables on which the biasing force is acting, has to represent every slow (i.e. associated to some metastability) degree of freedom in the system, in order for these methods to give meaningful free energy estimates [27]. Enhanced dynamics methods are somehow different from enhanced sampling methods such as umbrella sampling [28] and similar, which aims at performing equilibrium sampling in selected zones of the phase space using a time independent bias. The presence of a time dependent bias in accelerated dynamics methods introduces indeed some problems, as the results might depend on the bias build-up protocol. Recently, combinations of accelerated dynamics and accelerated sampling methods have been employed to get the best of the two approaches [10,29,30]. In these combined approaches, the phase space is first explored employing a dynamical biasing potential. After ergodicity is attained, build-up process is stopped, and the resulting (now, static) biasing potential is employed to perform equilibrium sampling in the way of umbrella sampling. The build-up scheme is arbitrary, but some algorithms, such as the combined metadynamics/umbrella sampling [29] and LEUS [10] approaches, uses a biasing potential updated by adding Gaussian terms at given intervals. The Gaussians are located at the instantaneous value of the CV in the metadynamics/umbrella sampling, while in the LEUS approach the CV space is discretized, and the Gaussians are located at the nearest grid element.

The free energy profile of the collective variables s is equivalently computed as

F (s) = -U bias (s) -k b T ln ρ bias (s), (4)
or by performing the re-weighting

exp[-F (s)/k b T ] = δ(s -s ′ ) exp[U bias (s)] b exp[U bias (s)] b . (5)
Here U bias is the biasing potential k b is Boltzmann constant, T is the temperature, ρ bias (s) is the probability density of finding the system at s during the umbrella sampling phase and . . . b denotes an average performed over the biased ensemble.

Simulation details

In this work we present results from different simulations (molecular dynamics, combined metadynamics/umbrella sampling in presence and absence of solvent, using CP or SP coordinates) involving β-d-glucose. In all simulations glucose molecules have been modeled using the gromos 45a4 force field [31], while the SPC model [32] has been employed for water molecules. Even though newer versions of the gromos force field provide a better description of the puckering properties of sugars [30,33], the 45a4 force field has been chosen for the present work because it allowed to estimate by standard molecular dynamics the equilibrium population of the chair conformers of glucose, therefore enabling us to check in a direct way the predictions of the coupled metadyamics/umbrella sampling approach with CP collective variables. All simulations has been performed at a constant temperature of 300 K, employing a Nosè-Hoover thermostat with time constant of 0.1 ps for the simulation in solution and a Langevin thermostat with time constant of 0.1 ps for the simulation in vacuum. In addition, for simulations in solution, the edges of the simulation box have been rescaled by means of the Parrinello-Rahman algorithm [34] with a time constant of 1 ps to keep a constant pressure of 1 atm. In this case, long range correction for dispersion forces have been applied. Bond lengths have been always kept constant at their equilibrium value by means of the shake algorithm [35], us- ing a relative tolerance of 0.1%. The simulation performed will be conventionally called simulation I-III. For simulations II and III, cutoffs of 1 nm have been applied to Lennard-Jones and electrostatic interactions, while for simulation I the cutoff length was 1.3 nm. In addition, electrostatic interactions have been screened by using a dielectric constant of 80, for the simulations in vacuum. Simulation I is a 150 ns long molecular dynamics simulation (with an integration timestep of 2 fs) of 512 glucose molecules in 25512 solvent molecules of solvent. This represent a big enough statistical sample to observe conformational changes, allowing to check the prediction of the metadynamics/umbrella sampling approach by estimating equilibrium population of chair conformers and, in addition, transition rates. Configurations along the trajectory were sampled every 10 ps and saved to disk for off-line analysis.

Simulation II is a metadynamics of a single glucose molecule in vacuum, performed by applying the bias potential on the two angular CP coordinates. During the metadynamics run, about 10 5 Gaussians were placed each every 0.1 ps, with a height of 0.12 kJ/mol, and a width of 0.05 rad in each angular direction. The free energy profile has been generated by calculating the value of the biasing potential on a 60 × 60 grid in the (θ, φ) space.

Simulation III is a metadynamics of a single glucose molecule in vacuum, performed by applying the bias potential on the three SP dihedrals. The number, height and frequency of Gaussians placed are the same as in simulation II, while their width was set to 0.2 rad in each of the three angular direction. The free energy profile has been generated by calculating the value of the biasing potential on a 60 × 60 × 60 grid in the (α 1 , α 2 , α 3 ) space.

CP and SP variables as reaction coordinates

The state-counting problem

Even though CP and SP coordinates are equally fit to describe puckered conformers from a purely geometrical point of view, differences can arise due to the shape of the puckering free energy hypersurfaces. Usually, in six-membered rings the CP puckering amplitude Q is not a slow degree of freedom, as it does not present metastabilities, and its free energy profile is such that conformers are generally populating only a thin spherical shell of the CP conformational globe [9]. It is therefore sufficient to introduce a biasing potential along the angular CP variables only, in order to get a proper exploration of the conformational phase space. Moreover, every coordinate system has to be complemented by a suitable assignment scheme for thermodynamic state, as geometrically different (in the sense of CP or SP coordinates) conformers located in the same free energy basin represent a single conformer in the thermodynamical sense.

Given the presence of only two variables, in case of CP variables visual inspection is sufficient to assign thermodynamic states, while quantitative schemes are required for SP variables. This can be easily seen by looking at the free energy profiles obtained from simulations II and III, presented in Fig. 2 and3, respectively. For simulation III, which employed SP dihedral angles, only some relevant cross sections of the free energy volume are presented, namely those at α 3 = -60, -35, 0, 35 and 60 deg. By looking at the free energy profiles presented, one can readily understand how different can be the task of identifying conformers location, free energy basins and barriers in the two coordinate systems.

In the LEUS investigation of Hünenberger and Hansen [10], the authors proposed three quantitative assignment schemes suitable for SP variables. Two of Free energy profile obtained from Simulation II. The superimposed solid lines defines the regions determined by the CD (black lines), projected onto the (θ, φ) plane . The conventional symbols of conformers associated to the different regions of the CD scheme are also shown. Blue and violet solid lines represent the projections of the ellipsoids for the MS scheme using a 5 and 10 kJ/mol threshold, respectively. these schemes, namely the closest distance (CD) and angle deviation (AD) ones are solely geometrical, while the minimum shape one (MS) employs a mixed energetic / geometrical criterion.

In the CD scheme, a conformation identified by the triplet of SP dihedrals β = (β 1 , β 2 , β 3 ) is assigned to the i-th state α i = (α i 1 , α i 2 , α i 3 ) if the quantity d i (β) = 1 -α i • β/( α i β ) is smaller than that of other states d j (β) (See Tab.1 for the definition of the triplets defining different puckered states). At a difference with the CD partitioning scheme, in the MS one states are identified by ellipsoids enclosing a local minimum and all the nearby conformations within an energy interval of 5 kJ/mol (the details are reported in [10]). Although the criterion is still partially geometrical, since it assumes an ellipsoidal shape of the basins, it makes use of the information from the free energy landscape, and is thus, among the three SP partitioning schemes, the closest to the visual inspection approach used with CP variables.

In the following we compare how three different schemes (visual inspection for CP coordinates; CD and MS for SP coordinates) partition the respective conformational space. It is therefore necessary to represent the thermodynamic states in a common framework. The natural way to do this is to project the volumes of states assigned in the SP space onto the CP plane. In Fig. 2 it is possible to see the different sub-volumes of the SP space corresponding to the 14 states determined by the CD assignment, projected onto the (θ, φ) plane. The 14 regions have been determined by reconstructing [6] the atomic positions of an ideal conformer (i.e. using bond lengths and angles of an ideal cyclohexane molecule). Starting from the CP puckering angles taken from a 500×500 grid on the (θ, φ) plane, SP dihedrals have eventually been calculated from the atomic coordinates, and the CD scheme has been applied (the computer code which computes atomic positions from CP parameters is made available at http://www.science.unitn.it/~sega/sugars.html.) In the background, the free energy profile of hydrated β-d-glucose obtained with the metadynamics/umbrella sampling approach [30] is shown. It is clearly seen that along the pseudorotational (equatorial) path, states determined by their thermally accessible conformational volume encompass more than one state defined by geometry considerations. This is the cases of the free energy basin around φ = 0 • that covers a group of five geometrically different conformers (namely B 25 , O S 2 , 3O B, 3 S 1 , B 14 ), or also the basin around φ = 180 • that groups four different conformers (namely 2 S O , B 3O , 1 S 3 , 14 B). This confirms what Hünenberger and Hansen noted in their work, namely, that purely geometrical criteria such as CD and AD do not allow to identify thermodynamic states (intended as free energy basins). The regions associated with the two chair conformers, at the north and south pole, respectively, seem on the contrary to capture quite nicely the features of the underlying free energy profile. Using the same procedure employed for the CD scheme, we calculated the projection of the MS ellipsoids obtained from the SP free energy profile onto the (θ, φ) plane, using two different threshold values of 5 and 10 kJ/mol, respectively. As it can be seen from Fig. 2, the MS scheme is able to describe noticeably well the shape of the main four free energy basins, however the full extent of the basins is not covered, at least with the 5 kJ/mol threshold. Notice that changing the threshold from 5 to 10 kJ/mol affects only the assignment of chair conformer basins, but not the boat and skew ones, as only the latter basins are less deep than 5 kJ/mol.

To check quantitatively the CD and MS schemes for chairs, we compared the percentage of these conformers along the evolution of simulation I estimated using the CD, MS and visual inspection schemes. Regarding the visual inspection, in particular, from the free energy profile reported in Fig. 2 one can see that the location of the free energy barrier for the chair is very well approximated by the region θ < 60. We reported in Fig. 4, the time evolution of the number of chair conformer, computed using the two criteria. The time evolution of population of chairs and inverted chairs (θ > 120) in the simulation box is practically indistinguishable for the CD and visual inspection schemes, and amounts to a relative root mean square difference between the two methods of 0.02 and 0.03% for chairs and inverted chairs, respectively, proving SP dihedrals and the CD approach quite robust for the identification of chairs. The MS scheme shows an interesting behavior: when using the 5 kJ/mol threshold, only roughly 60% of chairs (and inverted chairs) with respect to the visual inspection or CD schemes are observed, meaning a miscounting of states. However, since the measured populations are reduced by roughly the same ratio for both chairs and inverted chairs, the calculation of free energy difference using Boltzmann inversion formula gives a value of 10.0 kJ/mol, compatible with the other schemes (9.2 kJ/mol). This situation can be improved, as seen in Fig. 4, by using a 10 kJ/mol threshold, which enlarges the assigned basins of chairs. It has to be noted, however, that since the smaller basins of the boat and skew conformers do not change shape, the evaluation of their free energy difference will be indeed affected by the choice of the threshold.

From the results presented in Fig. 4 one can extract some additional, interesting information. A fit to an exponential decay, which assumes a reaction kinetic of the first order, shows that the population of chairs reaches the equilibrium value of 97.7% with a rate of about 13 ns -1 , showing still large oscillations on the 40 ns scale. The estimated population corresponds to a free energy of 9.3 kJ/mol, in good agreement with the estimate of 10 kJ/mol previously obtained by means of combined metadynamics/umbrella sampling [30].

Influence of the collective variable choice on convergence

Apart from the practical difficulties in identifying conformers discussed so far, one could argue that much more time should be spent to fill the free energy profile with the biasing potential in the SP case, with respect to the CP one, for the simple reason that the SP space is three-dimensional, whereas the CP is only twodimensional. In general, indeed, the time spent to sample a space grows exponentially with the number of its dimensions. However, this is true only if the different subspaces have the same characteristics, from an energetic and geometrical point of view. Including, for example, the total puckering amplitude in a metadynamics investigation with CP variables, will not increase dramatically the amount of phase space accessible at a certain energy, which will still have the shape of a spherical shell. Simple dimensional arguments are therefore not sufficient to estimate how the efficiency in filling the free energy profile depends on the set of collective variables in use. One possible estimate can be obtained, for example, by checking how many depositions are needed in order to fill with Gaussians a given amount of the energetically accessible portion of the free energy profile, ∆G, once the system has become ergodic.

Obviously, the portion filled by Gaussians for a given number of deposition will depend on the Gaussian height, but also on its width. Rescaling the number of depositions by the thermal energy in units of the Gaussians height is the obvious way to compare the outcome of two simulations that makes use of different heights. However, different coordinates requires in general the use of different values for the Gaussians width, but at a difference with the heights, they are in general not commensurable. In practice, a characteristic size always emerges naturally from the problem, and is the size of the smallest feature of interest in the free energy landscape. When setting a Gaussian width much smaller than this characteristic size, time will be spent without gathering important information on the free energy landscape, whereas if the width is too large, a metadynamics run will not be able to resolve properly the desired feature. In order to provide a meaningful comparison, the number of deposition have to be rescaled by the volume of a reference feature, in units of the width of the Gaussians. For the puckering free energy landscape of glucose one could choose the volume accessible within one unit of thermal energy, k B T from the minimum of a given conformer like, for example, 3O B.

In Fig. 5 we present the equivalent number of depositions (that is, normalized according to the procedure described above) needed to add a minimum number of points to each of the 14 states listed in Tab.1, identified using the CD scheme. The surface and volume of the 3O B conformer basin are 0.4 rad 2 and 1.2 × 10 -3 rad 3 , respectively. The tendency for both CP and SP cases is to reach a plateau value, which denotes an homogeneous filling of the free energy landscape. A fit to an exponential decay has also been performed, which gives an estimate plateau value of about 120 and 2300 equivalent Gaussians for the case of CP and SP variables, respectively. This result shows that using CP variables can lead to an tenfold increase in efficiency in filling the free energy landscape with respect to the SP dihedral angles, when using the same Gaussian height and an equivalent Gaussian width.

Conclusions

Different equivalent descriptions, such as the Cremer-Pople and Strauss-Pickett ones, can be adopted to characterize puckered ring conformers. Their equivalence on the pure geometrical level does not correspond, however, to an equivalence on the practical level, when the puckering coordinate are employed as collective variables in accelerated dynamics methods to compute puckering free energy landscapes.

The need for only two (angular) puckering coordinates to sample in an ergodic way the conformational space in the Cremer-Pople approach is at the root of two main advantages with respect to the usage of Strauss-Pickett dihedrals. The first advantage is that energetically-separated conformational basins are readily identified in the CP scheme, whereas when employing the three SP dihedrals, one has to rely on some algorithmic partitioning of the conformational space. Two of the checked partitioning schemes, namely the closest distance and the minimum shape, are shown to be prone to a misleading interpretation of the thermally accessible conformers of sugar rings. The reason for the miscounting in the closest distance scheme is related to its purely geometrical nature, while in the minimum shape scheme the miscounting is related to the choice of the energy threshold. The second advantage is that employing CP coordinates rather than SP dihedrals leads to a tenfold decrease of the time needed to homogeneously fill the conformational space.

Figure 2 .

 2 Figure2. Free energy profile obtained from Simulation II. The superimposed solid lines defines the regions determined by the CD (black lines), projected onto the (θ, φ) plane . The conventional symbols of conformers associated to the different regions of the CD scheme are also shown. Blue and violet solid lines represent the projections of the ellipsoids for the MS scheme using a 5 and 10 kJ/mol threshold, respectively.

Figure 3 .

 3 Figure 3. Cross sections (at α 3 = -60, -35, 0, 35and 60 deg) of the puckering free energy profile obtained from Simulation III.

Figure 4 .

 4 Figure 4. Time evolution of the population of 4 C 1 conformers in simulation I. Solid lines from top to bottom: (a) the CD scheme (almost coincident with the visual inspection scheme result); (b) the MS scheme with 10 kJ/mol threshold and (c) the MS scheme with 5 kJ/mol threshold. A fit to an exponential decay (dot-dashed line) for (a) is also presented.

  Gaussians / 10 3 minimum number of visits/10 3

Figure 5 .

 5 Figure 5. Efficiency of the CP and SP variables in filling the conformational free energy profile. Circles: CD scheme for SP dihedrals; Triangles: visual inspection criterion fo CP variables. The result of a fit to an exponential decay is also provided for each dataset.

Table 1 .

 1 Location of carbohydrate conformations in the SP dihedral space a .

	Conformer	α 1 (deg)	α 2 (deg)	α 3 (deg)
	4 C 1	-35	-35	-35
	1 C 4	35	35	35
	B 3O	30	-60	30
	B 25	30	30	-60
	B 14	-60	30	30
	3O B	-30	60	-30
	25 B	-30	-30	60
	14 B	60	-30	-30
	1 S 3	60	-60	0
	1 S 5	60	0	-60
	3 S 5	0	60	-60
	3 S 1	-60	60	0
	5 S 1	-60	0	60
	5 S 3	0	-60	60
	F o r
	P
	e
	e r
		R
		e
			v i e
				w
				O n l
				y

• ] and φ ∈ [0, 360 • ). The radial coordinate, Q, also called total puckering a From Rao et al.

[15]

.

September 3, 2010 20:9

URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics

Acknowledgements

The authors thank P. Hünenberger and H. Hansen for fruitful discussions. Use of the HPC cluster Wiglaf at the Department of Physics of the University of Trento is acknowledged. The authors are member of the Interdisciplinary Laboratory of


Jiří Škvor

Ivo Nezbeda

Percolation Line and Response Functions in Simple Supercritical Fluids

Keywords: Percolation line, Response functions, Supercritical fluids Percolation line, Response functions, Supercritical fluids

The question of a physical relevance (meaning) of the percolation line in simple supercritical fluids has been addressed using both extensive Monte Carlo simulations and accurate analytic equations of state. Thermodynamic and structural properties of two qualitatively different fluids, the Lennard-Jonesium (continuous model) and square-well fluid (stepwise model), have been studied in vicinity of their percolation lines over a range of pressures in order to find potential changes which may take place when an infinite cluster is detected in the system. Two different criteria for the occurrence of an infinite cluster, physical and configurational criteria, have also been used to assess their effect on the observed properties. It is found that the lines of extremes of various response functions run close to the percolation lines. Accounting for uncertainties in the definition of bonds and determination of the percolation line one may conjecture that extremes in some response functions occur when crossing the percolation line.

At certain thermodynamic conditions the molecules in bulk fluids may form aggregates (clusters) and/or a network whose occurrence affect considerably the behaviour of fluids. Formation or disruption of an infinite cluster of molecules (percolation transition) is then marked by a dramatic changes of certain properties of the system. Typical examples are a rapid change of electric conductivity of aqueous solutions [1,2], metal-nonmetal transition in supercritical fluids [3], or a sudden change of electric and dielectric properties of composites [4]. Percolation seems to be also a necessary condition for physical gelation [5].

Although the above examples refer to fluids, the concept of percolation originated in studies of lattice systems mimicking the crystalline phase for which the notions as, e.g., bond or the nearest neighbor are rather natural and for which also a number of theoretical results have been obtained [START_REF] Stauffer | Introduction to Percolation Theory[END_REF][START_REF] Sahimi | Applications of Percolation Theory[END_REF]. One of the most important of them is a universality which says that although the percolation threshold is system dependent, there are scaling functions and parameters which are universal for all systems of a given dimensionality and symmetry. However, when dealing with fluids (continuous systems referred to also as correlated systems) the situation is much less clear and a number of fundamental problems as, e.g., the definition of cluster itself, are encountered. Consequently, practically no theoretical results are available for fluids and various conjectures on percolation in fluids represent a mere extension of the results obtained originally for lattice systems. When dealing with fluids, the first problem to be addressed is the definition of cluster itself (i.e., a group of particles that are connected by an unbroken path of interparticle bonds) because there is no unambiguous criterion of particle connectivity. Then, in addition to determining the occurrence of clusters of various size, the ultimate goal of the percolation theory is to find the percolation threshold, i.e., the thermodynamic state at which a percolating (infinite) cluster occurs. Since molecular simulations are the only currently available tool the problem arises how to estimate the percolation threshold from simulations on finite systems. As regards the former problem, the most common way is to use a local concept, i.e., to define a cluster by means of bonds between pairs of molecules and the entire problem is then reduced to the problem of an appropriate definition of bonds. With the latter problem we dealt recently and concluded that, for supercritical fluids, the so-called wrapping probability (the cluster wraps the system if it is possible to get, starting from any particle of the cluster and moving along interparticle bonds, to an image of that particle in another replica [START_REF] Seaton | [END_REF]) provided a more accurate estimate of the percolation threshold when compared to the crossing probability (the cluster crosses the system if the maximal distance between some pairs of its particles is equal to or greater than the box length) [9,10]. Finally, we have also addressed recently the problem of potential universality for supercritical fluids and showed that the results exhibited not only strong temperature dependence but also dependence on the nature of the systems considered [11]. In other words, the percolation threshold in real fluids cannot be characterized by universal exponents of usual random percolating systems.

There are two other important issues related to supercritical fluids. The first one is the question where the percolation line terminates with lowering temperatures. Theoretical arguments suggest that it should terminate at the critical point [12] but vast majority of simulation results have not indicated any such tendency. Surprisingly, this important issue has in fact gained only little attention and we are aware of only two papers confirming the theoretical arguments for continuous models, for supercritical Lennard-Jonesium (LJ) [14,15]. With respect to the above mentioned ambiguities one reason may be an inappropriate choice for the bond definition or/and estimation of the percolation threshold.

There are two other important issues related to supercritical fluids. The first one is the question where the percolation line (i.e., the loci of the state points where the percolation transition/threshold occurs) terminates with lowering temperatures; in other words, what is the relation between the percolation and the critical point. This topic has been comprehensively reviewed by Coniglio [12]. With respect to the mentioned ambiguous criterion of the interparticle connectivity, there cannot be an unambiguous percolation line and, consequently, no its unique relation to the critical point. For example, it has been already well-known that application of a simple geometric definition of cluster tends to produce a percolation line ending below the critical temperature in the low-density region (e.g., [13]). On the other hand, we are aware of two papers which, using a 'proper' definition of cluster for a continuous model, supercritical Lennard-Jonesium (LJ), confirm that

the percolation line terminates at the critical point [14,15]. It is also known that the percolation line can be moved to the left or right side from the critical point by using different (more sophisticated) criteria (see, e.g., [16,17]). On top of that, another problem in interpreting the percolation line and its termination arises when an inappropriate method for the estimation of the percolation threshold is used.

The second issue is related to the fundamental question one should ask at the very beginning when studying percolation in supercritical fluids: Does the percolation threshold have any physical significance for supercritical fluids or is it only a pure mathematical construction extending otherwise meaningful results? Dramatic changes in the properties of fluids are observed primarily on mixtures (aqueous solutions, colloids) but what about pure supercritical fluids? This applies particularly to models with continuous intermolecular potentials. Based on physical intuition, one can understand that for models with a stepwise interaction potential as, e.g., the square-well (SW) potential or primitive models of associating or polar fluids [18][19][20], the occurrence of an infinite cluster may give rise to some changes in the thermodynamic or structural properties of the system. But what about continuous models as, e.g., the LJ potential? To our best knowledge, we are not aware of any result reporting observable changes in their properties related to the percolation. However, this does not seem to be the case of water which in fact has attracted a good deal of attention. Quite recently Bernabei et al. [21], combining neutron diffraction experiments with molecular simulation and the questionable assumption of universality, determined the percolation line for supercritical water and reported certain changes in the 3-dimensional arrangements of water molecules when crossing the percolation line. Partay and Jedlovszky [22] simulated recently supercritical SPC/E water and suggested that the percolation line might coincide with an 'imaginary' extension of the vapour-liquid (V-L) line in the pressure-temperature plane to the supercritical temperatures.

To summarize, the problem of the physical meaning of the percolation threshold in supercritical fluids, including water, has not been settled yet and the present paper attempts to address this issue and related problems. Percolation transition is not a genuine transition in the thermodynamic sense but may be viewed as a sort of a geometrical transition since a certain restructuring (rearrangement) does take a place when crossing the percolation line. We will therefore try to detect and characterize such changes and their relation to the thermodynamic behaviour of the fluid. Since the commonly used correlation functions are rather insensitive to tiny structural changes, we use the Voronoi polyhedra tessellation. As for thermodynamic properties, it is known that the so called thermodynamic response functions (derivatives of the state functions with respect to temperature, e.g., the isobaric heat capacity) [START_REF] Xu | Proc. Natl. Acad. Sci. USA[END_REF] exhibit extremes at a certain temperature, T m (P ), for pressures greater than the critical one, P c . This behaviour must/should be related to some structural changes in the system, potentially to percolation as suggested recently by Partay and Jedlovszky [22]. We will therefore focus on the response functions and extend their notion by considering the first-order derivatives also of the state variables (as, e.g., the isothermal compressibility) which exhibit extremes in temperature at a given pressure as well. After providing all the necessary definitions in the next section, we then determine the percolation line in two most common but qualitatively different simple fluid models, the LJ and the SW fluid. Further, two different definitions are used to consider two molecules to be bonded or not, the usual configuration criterion (also referred to as the Stillinger criterion) and the physical criterion (also referred to as the Hill criterion). Based on our own previous results [9,10] rence of an infinite cluster in the system and no assumption of universality is used. The obtained results are then summarized and discussed in section 4.

Basic definitions and computational details

The two common simple fluid models considered, the LJ and SW models, are defined by the following intermolecular potentials,

u LJ (r) = 4 σ r 12 - σ r 6 (1)
and

u SW (r) = +∞ for r < σ = - for σ ≤ r < λσ = 0 for r ≥ λσ , (2)
respectively. The depth of the potential, , and parameter σ are used henceforth to define the dimensionless temperature, T * = k B T / , pressure, P * = P σ 3 / , and density, ρ * = N σ 3 /V , where k B is the Boltzmann's constant, T is the temperature, P is the pressure, N is the number of particles, V is the volume of the system, and ρ is the number density. The range of the SW interaction has been set to λ = 1.5 which also defines the range within which two particles are considered bonded (configuration criterion). In the case of the LJ fluid, we choose from various possibilities for the configuration criterion for the limit of bonding the point of inflection of u(r), r b /σ = 26 7 1/6 , used already in previous studies, see e.g. [START_REF] Yoshii | [END_REF].

The energy-based Hill clusters are defined by means of the relative kinetic energy with respect to their interaction energy [25,26]: If v i and v j are the velocity vectors of particles i and j, resp., and m is their mass, then they are considered to be bonded if

m 4 (v j -v i ) 2 ≤ -u ij . ( 3
)
We used the common Metropolis Monte Carlo (MC) simulation in a N-V-T ensemble [START_REF] Allen | Computer Simulation of Liquids[END_REF] both for the location of the percolation threshold and for the investigation of other structural properties; for the application of the energy criterion the velocity components have been assigned to every particle at random from the Gaussian distribution characterized by the temperature of the system (for further details see the previous paper [11]). The correlation functions and the structure factors were determined at several state points around every detected percolation threshold point and further structural analysis was carried out by means of the Voronoi tessellation [START_REF] Medvedev | [END_REF]29] over a wide range of densities. The extension of the V-L line to the supercritical temperatures has been determined by an extrapolation of the Clausius-Clapeyron equation ln P = a -b/T , where a and b are fitted parameters for the model considered.

The response functions we consider here are the coefficient of thermal expansivity,

α P = 1 V ∂V ∂T P = - 1 ρ ∂ρ ∂T P , (4)
the coefficient of isothermal compressibility,

κ T = - 1 V ∂V ∂P T = 1 ρ ∂ρ ∂P T , (5)
and the isobaric heat capacity,

c P = 1 N ∂H ∂T P , ( 6
)
where H is enthalpy. These functions exhibit various extremes or even singularities (as, e.g., the heat capacity at the critical point) or change sign (as, e.g., the thermal expansivity which indicates an anomaly in the isobaric temperature dependence of density). The loci of the extremes of the response functions as functions of temperature at constant pressures form in the P -T plane a line which we are going to refer to as the 'response function extremes line (RFEL)'. To determine the response functions and RFEL's, for the LJ fluid we use the most accurate equation of state (EOS) available to date, Kolafa-Nezbeda equation [30], and for the SW fluid a generalized EOS developed by Patel et al. [31].

Results and discussion

Following the strategy outlined above we determined from simulations two types of the percolation line corresponding to the two different definitions of bonds (configurational and physical) and the results are shown in figures 1 and 2. As it is seen, both lines follow a different course. Although it may seem that at least the Hill's percolation line tend to go to the critical point, the detailed examination shows that this is not the case. And as for the Stillinger's percolation line, it evidently goes away from the critical point. As already mentioned, the primary property of the fluids from onto which all other properties reel off are superimposed is the structure. As regards structural characteristics, we have not found any significant (reasonably detectable) changes in the course of either correlation functions or the structure factors that might indicate a structural rearrangement. From a number of various characteristics of the Voronoi tessellation (volumes of Voronoi polyhedra, etc.) we examined, the most interesting seems the skewness of the distribution of the number of nearest neighbors (in terms of the Voronoi tessellation) for one simple reason: this function exhibits, as a function of density ρ at a given temperature T , a maximum, see figure 3, and the line connecting these extremes (skewness extreme line; SEL) originates from (vicinity of) the critical point.

Examination of figures 1 and 2 reveals two different patterns of the behaviour not in the dependence on the fluid model but in dependence on the thermodynamic property considered. Quantities directly related to the structure are the coefficient of isothermal compressibility and the skewness of the distribution of the number of nearest neighbors. As we see, the line of extremes of the former quantity runs very close to the SEL and both also follow (approximately) the course of the percolation line based on the configurational definition all along the considered temperature range. The former coincidence expresses the fact that the coefficient of isothermal compressibility reflects the response of density to the applied pressure which may give rise, at certain conditions, to a structural rearrangement (repacking) of the molecules and the same applies to the skewness. Since the percolation line involves explicitly the bonding in the geometrical sense (proximity of molecules), it is understandable that the structural rearrangement may simultaneously give rise to establishing of an infinite network, i.e., percolation. Quite different pattern of behaviour we find in the case of the other two RFELs. First, the RFEL's of the isobaric heat capacity and of the coefficient of thermal expansivity are not monotonous functions of pressure and may thus hardly be related to any percolation line, which are monotonous, over a wide range of pressures/temperatures. Second, both these functions are derivatives with respect to temperature and involve thus, implicitly, the effect of the change in energy. If any percolation line is to be related to them then it is the percolation line based on the Hill's criterion which is thus in this case more appropriate. In fact, for the SW fluid this percolation line runs, in the vicinity of the critical point, quite close to their RFELs.

The last curve shown in figures 1 and 2 is the extension of the V-L line. In principle, this line does not have any physical meaning and, as we see, it cannot be related to any pattern of the behaviour of the quantities shown in these figures.

As a counterpart to the P -T diagram used in discussions of any transition, we show in figure 4 the above discussed quantities also in the T -ρ plane which is a common practice when dealing with percolation. In this case it is even more difficult to find a common pattern for the percolation lines and other characteristics of the considered fluids. In fact, the results are different for each fluid with the exception of the percolation line given by the configurational definition of the bond: In both cases the percolation line ends anywhere but in the critical point. Some relation between the Hill's percolation line and the energy dependent quantities (α P and c P ) can be detected for the SW fluid over an extended range of the supercritical temperatures but no such pattern is observed for the LJ fluid.

Conclusions

Percolation is a phenomenon which is known to affect considerably the properties of the solid phase but whether the same applies also to the fluid phase remains still an open question. We have investigated the behaviour of certain thermodynamic and structural properties of two simple supercritical fluids and found out two different types of their behaviour, both with respect to the lines of percolation and among themselves.

Arguing that the percolation in the fluid is connected, in some way, to a structural rearrangement of the molecules in the system, then the configuration criterion seems the most appropriate. The thermodynamic property reflecting such a structural change is the coefficient of isothermal compressibility whose RFEL (the line connecting the state points at which the given response function, as a function of temperature T at constant pressure P , reaches its maximum) follows very closely the course of the SEL (the line connecting the state points at which the skewness of the distribution of the number of nearest neighbors in terms of the Voronoi tessellation, as a function of density ρ at constant T , reaches its maximum) and both follow, approximately, also the course of the configuration percolation line in the P -T plane. The other case involves properties which are, either explicitly or implicitly, associated with changes in energy, the coefficient of the thermal expansivity and the isobaric heat capacity. Their RFELs fall into the region of the percolation line defined by the Hill's criterion. However, with respect to the behaviour of the RFEL's no simple conclusion can be drawn.

Some authors [22,[START_REF] Xu | Proc. Natl. Acad. Sci. USA[END_REF] in supercritical fluids may be an extension of the vapour-liquid (V-L) line. However, we have not found any evidence for this surmise, except for very close proximity of the critical point; in fact, outside the critical region the V-L line runs away from both the percolation lines and also away from the RFELs and SEL. Nonetheless, with respect to ambiguity of the definition of the percolation line [17] such a connection cannot be, in principle, excluded.

To summarize, we have attempted to find a link between the percolation and observable changes in the thermodynamic behaviour of supercritical fluids considering two types of simple fluids. With respect to the mentioned ambiguity in the definition of cluster (and percolation) in fluids, only indications of such a potential link between the percolation and extremes of (at least some) response functions have been found, at least in the vicinity of the critical point. For a more definite result at least other classes of fluids must be also considered, particularly associating fluids (e.g., water and methanol) and, for completeness, also strongly polar fluids (e.g., aceton). Thus, the problem of the physical relevance of the percolation line in supercritical fluids remains still open and deserves further investigations. A revised version of the above specified paper has just been uploaded. The referee asked for a minor revision and the paper has been revised accordingly.

Most changes take place in the Introduction. Specifically: 1. Wrapping and crossing probabilities are explicitly defined. 2. Percolation line is defined.

3. As for the term 'Response function', which is clearly defined in the paper, this is not our `invention'. It's used in the field and a reference has been added for further details. 4. One paragraph dealing with the connection between the percolation line and the critical point has been completely reformulated and extended, and additional references added.

Structural properties are discussed in terms of the skewness for reasons which are discussed in the paper in details.

The acronyms used in the Conclusions has been detailed again for better readibility.

All the changes to the MS are highlighted by color: the green color indicates the new text, red color the text which must be deleted.

I believe that after these changes the paper will be acceptable for publishing in Mol Phys.

With best regards,

Figure 1 .

 1 Figure 1. Percolation lines, loci of the extremes of the response functions, and the skewness extreme line in the supercritical Lennard-Jones fluid.

Figure 2 .

 2 Figure 2. The same as figure 1 for the square-well fluid.

Figure 3 .

 3 Figure 3. The skewness, γ 1 , of the distribution of the number of nearest neighbors (in terms of the Voronoi tessellation) as a function of density at given temperatures for the Lennard-Jones (left) and square-well (right) fluids.

Figure 4 .

 4 Figure 4. Percolation lines, loci of the extremes of the response functions, and the skewness extreme line in the supercritical Lennard-Jones (left) and square-well (right) fluid.

Page 11 of 11 URL

 11 : http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics

  , we consider only wrapping clusters to identify the occur-

	eptember 2, 2010	20:41	Molecular Physics	SkvorNezbeda˙Percolation
		4		Jiří Škvor and Ivo Nezbeda
			F o r
			P
				e
				e r
				R
				e
				v i e
				w
				O n l
				y

  argue that the line which may indicate a structural change

	eptember 2, 2010	20:41	Molecular Physics	SkvorNezbeda˙Percolation
				Molecular Physics	7
			F o r
			P
				e
				e r
				R
				e
				v i e
				w
				O n l
				y

Acknowledgement

This work was supported by the Grant Agency of the Academy of Sciences of the Czech Republic (Grant No. IAA200760905).