Thomas Treal

Philip L Jackson

Aurore Meugnot
email: aurore.meugnot@universite-paris-saclay.fr

Combining trunk movement and facial expression enhances the perceived intensity and believability of an avatar's pain expression

Introduction

The last decade has witnessed a great enthusiasm for the use of computer-generated (CG) characters (avatars) in social and affective neuroscience studies. This type of technology promises to meet the idealistic compromise between controlled and ecological paradigms (Bombari, Schmid Mast, Canadas & Bachman, 2015;[START_REF] Pan | Why and how to use virtual reality to study human social interaction : The challenges of exploring a new research landscape[END_REF][START_REF] Wykowska | Embodied artificial agents for understanding human social cognition[END_REF]. That being said, the development of platforms allowing the animation of avatars remains complex and has given rise to the key issue of their realism (i.e., the degree of human likeness of the virtual character), which may be crystallized through the 'fear' of falling into the Uncanny Valley phenomenon [START_REF] Kätsyri | Virtual Faces Evoke Only a Weak Uncanny Valley Effect : An Empirical Investigation With Controlled Virtual Face Images[END_REF][START_REF] Macdorman | Too real for comfort? Uncanny responses to computer generated faces[END_REF][START_REF] Mathur | Navigating a social world with robot partners : A quantitative cartography of the Uncanny Valley[END_REF][START_REF] Reuten | Pupillary Responses to Robotic and Human Emotions : The Uncanny Valley and Media Equation Confirmed[END_REF]; see [START_REF] Brenton | The Uncanny Valley: does it exist and is it related to presence[END_REF][START_REF] Kätsyri | A review of empirical evidence on different uncanny valley hypotheses: support for perceptual mismatch as one road to the valley of eeriness[END_REF] for critical reviews). This well-known phenomenon refers to 'a profound sense of unease provoked by an almost but not fully realistic character that prevents any sense of presence or emotional engagement' (Mori, 1970(Mori, /2012). In the domain of emotion research, the use of anthropomorphic avatars may be crucial. If the avatar is perceived as a simple machine, its facial emotional expression would not be believable. It would be considered as fake and artificial and would unlikely trigger any spontaneous emotional reaction [START_REF] Bogdanovych | What makes virtual agents believable?[END_REF]. Such an experimental paradigm would not be useful in the context of inferring the mechanisms underlying others' emotion processing and regulation.

Scholars argued that the lack of believability towards the avatar's behavior, and the eerie feeling that may accompany it, would be elicited by a perceptual mismatch in the realism of the avatar's features, with some appearing human and others very artificial (e.g., artificial eyes on a realistic face or disparity between voice and face realism levels, [START_REF] Macdorman | Too real for comfort? Uncanny responses to computer generated faces[END_REF][START_REF] Kätsyri | A review of empirical evidence on different uncanny valley hypotheses: support for perceptual mismatch as one road to the valley of eeriness[END_REF]. In particular, realism inconsistency could result from a gap between a highly graphically realistic avatar and a low realism in terms of motion and behavior [START_REF] Brenton | The Uncanny Valley: does it exist and is it related to presence[END_REF]; de Gelder, [START_REF] De Gelder | Virtual reality and the new psychophysics[END_REF]. This discrepancy could be the result of the disproportional amount of efforts put into pictural (or graphical) realism in avatar development [START_REF] Brenton | The Uncanny Valley: does it exist and is it related to presence[END_REF] as opposed to their movements -perhaps due to the difficulty of gesture animation (Hartmann, Mancini & Pelachaud, 2005;Normoyle, Lui, Kapada, Badler & Jorg, 2013). This is nevertheless quite surprising considering the bulk of evidence supporting the role of body expression in conveying emotion [START_REF] De Gelder | Towards the neurobiology of emotional body language[END_REF][START_REF] De Gelder | The perception of emotion in body expressions : Emotional body perception[END_REF]. It is even more surprising when one considers that Mori -who is extensively cited in the field of affective computing -attributed the highest level of affinity with an avatar that portrays a moving healthy character (Mori, 1970;2012). One might expect this observation to instill momentum into the development of both the pictural and behavioral features of virtual characters, notably when using them as stimuli to explore the interactive mechanisms underlying emotion perception and regulation.

To our knowledge, studies exploring the influence of body motion on the perception of facial emotions 'expressed' by an avatar are scarce [START_REF] Buisine | The Role of Body Postures in the Recognition of Emotions in Contextually Rich Scenarios[END_REF][START_REF] Clavel | Combining Facial and Postural Expressions of Emotions in a Virtual Character[END_REF][START_REF] De Melo | Influence of autonomic signals on perception of emotions in embodied agents[END_REF][START_REF] Wu | Effects of virtual human animation on emotion contagion in simulated inter-personal experiences[END_REF]. For instance, in an elegant study, De [START_REF] De Melo | Influence of autonomic signals on perception of emotions in embodied agents[END_REF] demonstrated that implementing respiratory patterns corresponding to various emotions (excitement, pain, relief, anger, fear, panic, boredom and startle) to an avatar increased the perception of these emotions in a virtual character, compared with a control (avatar without body animation) condition. In another seminal work, Clavel and collaborators (2009) examined how body posture may contribute to the recognition of basic emotion (fear, anger, joy, surprise, sadness) expressed by a virtual character. The authors showed that participants recognized these emotions more easily when the avatar's face and posture were jointly presented, compared to conditions with either face or posture alone [START_REF] Clavel | Combining Facial and Postural Expressions of Emotions in a Virtual Character[END_REF]. The same group replicated similar results in scenarios of an avatar expressing emotions while telling a story including ironic and narrative items [START_REF] Buisine | The Role of Body Postures in the Recognition of Emotions in Contextually Rich Scenarios[END_REF]. The authors showed that the negative emotions were judged to be more intense when the avatar adopted a posture in congruence with its facial expression compared to a condition where the avatar was still. Interestingly, Buisine and collaborators (2014) also examined the influence of an idle posture (i.e. small variations in body postures or balance changing) on participants' perception of the avatar's emotions, compared with a fixed posture.

Idle motion was coined by [START_REF] Egges | Personalised real-time idle motion synthesis[END_REF] to describe neutral (i.e. non-communicative) movements, including "changing balance because of fatigue, small variation in body posture caused by small muscle contraction, eye blinking". In other terms, idle motion refers to subtle movements that are not always perceptible if one pays no attention to them. Such movements are almost always present in humans [START_REF] Perlin | Real time responsive animation with personality[END_REF]. These movements are generally lacking in avatar animation, but they might be a core ingredient to creating the illusion that the avatar is alive. There are some exceptions, for instance idle motion implemented by default in a conversational avatar in [START_REF] Groom | Evaluating the effects of behavioral realism in embodied agents[END_REF], and the control condition in De [START_REF] De Melo | Influence of autonomic signals on perception of emotions in embodied agents[END_REF]. [START_REF] Buisine | The Role of Body Postures in the Recognition of Emotions in Contextually Rich Scenarios[END_REF] observed that participants tend to recognize more emotions in the idle than in the still condition, while there was no significant difference between the two conditions for the perceived emotion intensity, as well as the evaluation of emotion realism.

Although mixed, these results are important to consider because, to our knowledge, this is the first study examining how this type of body motion -which may be a critical component of the avatar's animation -may influence the perception of emotion as well as its overall realism.

Moreover, other authors showed that participants' emotional reactions were stronger when viewing the pain of a virtual patient whose medical condition was growing worse in a dynamic condition (the avatar is moving) compared to a condition where the virtual patient remained static [START_REF] Wu | Effects of virtual human animation on emotion contagion in simulated inter-personal experiences[END_REF]. This latter study may be of particular interest for the present work, which has used a similar virtual reality platform using recently avatars expressing pain in a medical context [START_REF] Tessier | Toward dynamic pain expressions in avatars : Perceived realism and pain level of different action unit orders[END_REF]. Finally, authors investigated the psychophysiological response to emotional stimuli by the means of pupil size changes [START_REF] Bradley | The pupil as a measure of emotional arousal and autonomic activation[END_REF][START_REF] Kret | Emotional signals from faces, bodies and scenes influence observers' face expressions, fixations and pupil-size[END_REF][START_REF] Reuten | Pupillary Responses to Robotic and Human Emotions : The Uncanny Valley and Media Equation Confirmed[END_REF]. [START_REF] Bradley | The pupil as a measure of emotional arousal and autonomic activation[END_REF] observed that pupil diameter was larger when participants were presented with pictures eliciting both pleasant and unpleasant contents than pictures with a neutral content. According to the authors, pupil dilation may be caused by sympathetic nervous system activity and reflect emotional arousal, which is a dimension of emotional reaction. Indeed, emotional reaction is often defined as the combination of two systems: the emotional valence (pleasantness/positive vs unpleasantness/negative), which indicates the direction of the emotional reaction, and the emotional arousal, which refers to the intensity of the emotional reaction [START_REF] Lang | Emotion, attention, and the startle reflex[END_REF][START_REF] Russell | Affective space is bipolar[END_REF]. Thus, pupil dilation is well established evidenced as a continuous and objective marker of emotional arousal when viewing facial expressions. Other authors recently showed that pupillary responses may also change when observing a virtual character's facial expression . Pupil responded less vigorously if participants judged the avatar's facial expression as eerie or unattractive [START_REF] Reuten | Pupillary Responses to Robotic and Human Emotions : The Uncanny Valley and Media Equation Confirmed[END_REF].

Thus, two studies were conducted. Experiment 1 aimed to determine the impact of different levels of trunk movement on a) the perceived intensity of pain expressed by an avatar, b) the believability of its facial pain expression, and c) the participant's emotional arousal assessed by pupil size recording. It was hypothesized that the implementation of trunk movement with a low amplitude (i.e., idle motion representing postural oscillation, typical of human dynamic equilibrium) would improve participants' perception of the avatar's pain and believability compared to a static condition. Larger amplitudes of trunk oscillation (i.e. trunk rocking) were supposed to express a pain behavior and would strengthen the congruency between face and body. Consequently, the avatar's pain expression would be intensified, up to a certain point where trunk oscillation would seem exaggerated and might make the avatar's pain seem weird, so less believable. Moreover, the pupil would dilate more if the avatar's pain were perceived to be more believable and more intense, namely in the conditions where trunk movement has been added to it compared with the static condition.

Experiment 2 using psychophysical cues was conducted to objectively classify the different conditions of trunk movement as idle motion or trunk rocking. The results observed in the Experiment 1 regarding the impact of trunk movement on participants' perception of the avatar's pain (Pain Intensity and Pain Believability) and the emotional arousal elicited by pain (Pupil Size Change) will be further discussed in light of the categorization obtained in Experiment 2.

2.

Experiment 1: Effect of body movement added to an avatar on the way one perceives facial expression of pain.

Method

Participants

Twenty-four participants (mean age = 22.04 years, SD = 3.61, range = 18-33 years; 13 men (54.2%)) were recruited through advertisements sent via South-Paris University's student mailing lists or posted on the campus of South-Paris University. Exclusion criteria included any reported history of painful, neurological, psychiatric disorder, work experience as a healthcare provider and an uncorrected visual impairment. This study was approved by the Ethics Committee of Paris-Saclay University (#CER-Paris-Saclay-2018-16) and participants gave written informed consent to take part in this study.

Procedure

Participants sat comfortably on a chair with their head supported by a chin-rest and facing a 20-inch monitor (60Hz refresh state, resolution 1680 x 1050 pixels) at a viewing distance of 55cm. A 9-point calibration of the eye-tracker was conducted until a good calibration was accepted (i.e. with a maximum error angle of 1°). After the calibration phase, participants were instructed to not move their head until the end of the experiment. Then, they completed the Pain Intensity task and the Pain Believability task. The order of the tasks was counterbalanced between the participants. After the experiment, a debriefing was held to inform the participants of the true objective of the study.

Material and tasks

Seven video clips (duration = 4 seconds) displaying the upper body from the trunk up of a male avatar originated from the Empathy-Enhancing Virtual Evolving Environment (EEVEE) [START_REF] Jackson | EEVEE : The Empathy-Enhancing Virtual Evolving Environment[END_REF] were created with Blender® (BlenderFoundation, Amsterdam, Netherlands). The avatar facing the camera at an 86.8° angle showed a dynamic expression of pain, either with or without an oscillatory movement of the trunk.

The avatar's pain expression was created using the Action Units (AUs) according to the Facial Action Coding System (FACS) [START_REF] Ekman | Facial action coding systems[END_REF], which are typical of pain expression [START_REF] Prkachin | The consistency of facial expressions of pain: a comparison across modalities[END_REF]: AU4 (brow lower), AU6 (check raiser), AU7 (lid tightener), AU9 (nose wrinkle), AU10 (upper lip raiser) and AU43 (eyes closed); and five occasionally AU27 (mouth stretch). The intensity of the pain expression was fixed at 60% of the maximal contraction of the AUs described above. Each video began with the avatar showing a neutral expression (sets of AUs at 0). The level of AUs linearly increased to reach 60% of contraction for 1 second. The pain expression was maintained for 3 seconds until the end of the clip (Figure 1). In six video clips, the avatar's trunk was moving to and fro at different amplitudes (0.33°; 0.67°; 1°; 1.67; 2°; 3° conditions). Trunk movement was limited to the anterior-posterior part and the degree of oscillation was determined to be at some point between the center of the mass and the top of the head. In each video, the avatar first stood perpendicular to the ground and moved forward (half amplitude), then moved back (total amplitude) and finally returned to its initial position (half amplitude) (Figure 2). In the seventh clip , no trunk movements were added to the avatar (still condition). imaginable"]). Each clip was displayed 10 times during the experimental phase. In total, this phase was composed of 70 randomized trials. A trial began with a centered fixation cross which was displayed for 1500ms. Then, the video clip was presented (duration=4 seconds), followed by the VAS. Every 10 trials, a break-screen was displayed until participants pressed the space key to continue the task. Before the experimental phase, the participants performed a familiarization phase composed of three trials with a break-screen between the second and third familiarization trials.

The Pain Believability task. Participants had to judge whether the avatar displayed in the clips was genuinely experiencing pain or not. To this end, they were informed that the avatar was animated with a pain expression captured either on individuals who were actually experiencing pain ("True pain") or actors imitating a pain behavior ("Fake pain"). Like in the pain intensity task, this task was composed of a familiarization phase with three stimuli and a break-screen, followed by an experimental phase of the seven video clips presented 10 times at random (70 trials in total). A trial began with a centered cross fixation lasting 1500 ms.

Then, a video clip was presented followed by a screen displaying two panels: "True pain" and "Fake pain", respectively centered on the left and the right part of the screen [in French:

"VRAI" and "FAUX"]. To respond, participants had to click on one panel with a mouse.

Eye-tracking

Pupil size was collected continually during the experiment using a monocular desk-mounted eye-tracker (EyeLink 1000, SR Research, Canada) at a rate of 1000 Hz and registered in arbitrary units using Experiment Builder (SR Research, Canada).

Data Analysis

The Pain Intensity task. The position of the cursor along the scale was converted to numerical values between 0 ("no pain") and 100 ("worst pain imaginable"). The mean Pain Intensity was calculated for each condition of trunk movement (condition 0° included).

The Pain Believability task. The Pain Believability rating corresponding to the percentage of "True pain" responses was calculated for each condition of trunk movement (condition 0° included).

Eye-tracking. Pupil size was analyzed using DataViewer (SR Research, Canada). One participant was removed from the analysis because of a high data loss during eye tracking recording (37.5%, mean date loss for others: 1.7% (SD=1. 0.67°; 1°; 1.67; 2°; 3°), with alpha set at 0.5. Pairwise comparisons were applied between conditions with a significant difference from 0 using Bonferroni correction if required.

Results

The Moreover, there was no significant difference between conditions 1°, 1.67°, 2° and 3° (p > .378), indicating the reach of a ceiling.

Summary of Experiment 1

Participants rated the avatar's pain expression as gradually more intense (condition 0°: hypothesis, the largest amplitude of trunk oscillation (3°) did not provoke a decrease in the perceived intensity and believability of pain in the avatar compared with the still and lower conditions of trunk movement. Unexpected results were also observed for pupillary response . The pupil dilated more strongly only when the avatar oscillated at a large amplitude (conditions 1°: 1.71%, 95% CI [0.41, 3.01] and 3°: 3.28%, 95% CI [0.85, 5.71]) than when it stayed still. Surprisingly, there was a lack of significant change for the conditions 1.67° and 2° compared with the static condition, while there was a significant change for the preceding (1°) and next (3°) conditions (Figure 5). Importantly, the lower amplitude of trunk movement (condition 0.33°) did not differ from the still condition for the three dependent variables (i.e.

mean
Pain Intensity, Pain Believability and pupil size). This lack of difference may suggest that this amplitude was too low to be perceived, and raised the question of its definition as idle motion.

A second study (Experiment 2) was conducted to objectively classify the different conditions of trunk movement as idle motion or trunk rocking, which may reflect a pain behavior when combined with facial pain expression. The results observed in Experiment 1 regarding the impact of trunk movement on participants' perception of the avatar's pain and the emotional arousal elicited by pain will be further discussed in light of the categorization obtained in Experiment 2.

Experiment 2: Categorization of trunk movement as idle motion or trunk rocking

A classical psychophysical tool was used to determine the categorical boundaries between the different conditions of trunk movement [START_REF] Bonnet | Manuel pratique de psychophysique[END_REF][START_REF] Courbalay | Big Five personality traits contribute to prosocial responses to others' pain[END_REF][START_REF] Wichmann | The psychometric function: I. Fitting, sampling, and goodness of fit[END_REF]. The aim was to classify them as idle motion or trunk rocking, which may reflect a pain behavior in the context of this study (i.e. when combined with facial pain expression). The video clips used in Experiment 1 were presented to a group of naive participants in a visual discrimination task. Herein, the participants were asked whether or not they had perceived the avatar's chest moving in each video clip. To be considered an idle motion, a trunk movement must be perceived (i.e. above a detection threshold), yet not systematically (i.e., under a systematic detection threshold), as these movements are by definition subtle and thus difficult to detect. The boundaries of idle motion were then calculated by fitting a psychophysical function to the detection rating of each amplitude of trunk movement (see below section 3.1.2 for an explanation of the related psychophysical indices). The amplitudes of trunk oscillation which fell above the systematic detection threshold were considered as the bodily expression of a sustained pain.

Method

Participants

The participants were twenty young adults recruited in the Faculty of Sport Science of

Material and Task

The same video clips as in the first experiment were presented to participants sat in front of the computer screen (20-inch, 60Hz refresh state, resolution = 1680 x 1050 pixels) at a distance of 55 cm from the head. The clips were presented 10 times per condition (70 randomized trials). A trial began with the display of a clip, and after participants had to respond whether the avatar's trunk moved or not by clicking "YES" or "NO" [in French, "OUI" or "NON"] labels, respectively presented on the center-left and center-right of the screen .

A logistic function y = 1/(1+exp(-(x-PSE)/JND)) was fitted to the average of the percentage of "YES" responses for each condition. The Point of Subjective Equality (PSE) corresponds to the detection threshold (the lowest amplitude of trunk movement that participants are able to detect). The Just Noticeable Difference (JND) refers to the precision of categorization. It corresponds to the minimum amplitude necessary to perceptually differentiate two sizes of oscillations. The JND is also the amplitude that must be added to the PSE to elicit a quasisystematic detection (PSE + JND = quasi-systematic detection threshold).

Results

As illustrated in Figure 6, the PSE was associated with an amplitude of trunk oscillation of 0.57° and the JND was equal to an amplitude of 0.20°. Our results indicate that participants were able to detect oscillation from 0.57° of oscillation (detection threshold) and that an amplitude of 0.77° (=PSE + JND) of oscillation led to a quasi-systematic detection.

Summary of Experiment 2

The results show that trunk movement at an amplitude of 0.33° was not detected by participants (under the PSE=0.57°). Participants likely did not perceive the avatar's trunk oscillation in the condition 0.67°, which was situated under the systematic detection threshold (PSE + DJP = 0.77°). Moreover, the larger amplitudes of trunk oscillation (from 1° to 3°)

were systematically detected.

Thus, except for the condition 0.33° which cannot be classified as idle motion, these results confirmed the categorization that was suggested in Experiment 1, namely that the condition 0.67° may refer to idle motion and the conditions 1°, 1.66°, 2°, 3° may be interpreted as the body expression of a sustained pain when combined with facial pain expression.

General Discussion

In the field of affective and social neuroscience, the increasing use of naturalistic, yet well-controlled CG characters (avatars) has led to a new multidisciplinary research avenue working towards the development of anthropomorphic characters, which not only look fully human but also convey 'genuine' emotions. The present work tested whether body motion combined with facial expression increased the perceived intensity and the sense of believability towards the pain expressed by such CG characters, i.e. an avatar from a platform that is specifically dedicated to the study of empathy [START_REF] Jackson | EEVEE : The Empathy-Enhancing Virtual Evolving Environment[END_REF]. Moreover, pupil size was recorded during the interaction with the CG character to assess whether pain expression of the avatar elicited a greater emotional arousal when body animation was there.

Effect of trunk oscillation as idle motion on participants' perception of the avatar's pain

Results of Experiment 1 showed that the condition 0.33° did not modulate participants' perception of the avatar's pain, as indicated by the lack of significant difference between this condition and the static condition, either on perceived Pain Intensity and Pain Believability or on emotional arousal inferred by pupil size. Results of Experiment 2 confirmed that the oscillation at 0.33° was not detected by the participants and, thus, was too low to be considered as idle motion.

By contrast, the condition 0.67° was likely inferred as idle motion since it corresponded to its definition: the oscillation was detectable as it was largely above the detection threshold. However, it was too subtle to be systematically perceived by the participants, even when their attention was focused on the avatar's trunk movement.

Importantly, results for Experiment 1 revealed that the facial pain expression of a CG character was perceived to be more intense and more convincing (i.e. judged more frequently as actually having been experienced) in the presence of this slight trunk oscillation than when the avatar's body was still. It is likely that the avatar looked more 'alive' when it was animated with subtle trunk movement than when it remained still, providing more consistency between graphical and behavioral realism. In fact, the 'statue-like' behavior of the avatar would mismatch with the dynamic facial expression in the static condition. Thus, the facial pain expression would convey less emotion, resulting in a lower perception of pain intensity and believability than in the idle condition.

Reuten and collaborators (2018) have previously shown that the pupil responds less to emotional facial expressions of robots that were judged as uncanny by the participants, compared to more human-like robots. In concordance with this study, it was hypothesized that if the avatar's facial pain expression was judged to be more intense and believable in the idle condition, it would also evoke a greater pupil response. However, the change in the avatar's pain perception was not accompanied by a greater emotional arousal when low trunk movement was added, as indicated by a similar pupil size between the static (0°) and the idle (0.67°) conditions. The lack of idle effect on pupil size may result from the very simple, yet quite artificial method used to animate the avatar's trunk. Indeed, postural oscillations down to a to and fro movement in the anterior-posterior plan (see section 2.1.3.2. Material and tasks), while they referred to a much more complex pattern of motion, often spatially represented by an ellipse [START_REF] Schubert | Ellipse area calculations and their applicability in posturography[END_REF]. Moreover, although postural oscillations are generally associated with respiration, which directly influences postural control [START_REF] Hunter | Respiratory components of human postural sway[END_REF], breathing movements were not included here. Finally, the avatar's trunk movement was not natural, since it was applied 'mechanically' into the animation software. Motion capture would have been much more suitable as it would take into account natural and unconstrained movements, which are well-known to improve the avatar's acceptability [START_REF] Pan | Why and how to use virtual reality to study human social interaction : The challenges of exploring a new research landscape[END_REF][START_REF] Piwek | Empirical evaluation of the uncanny valley hypothesis fails to confirm the predicted effect of motion[END_REF][START_REF] Thompson | The Perception of Humanness from the Movements of Synthetic Agents[END_REF]. Moreover, a recent study demonstrated that natural (i.e., biological) motion produced greater pupil dilation than mechanical motion does (William, Cristino & Cross, 2019), which supports the use of natural motion when assessing emotional arousal elicited by a CG character. However, avatar movements with motion capture would need to be pre-recorded and relayed, which implies constraints to a potential use of avatars in an interactive and flexible manner.

Overall, even though the method used for body animation limits the ecological scope of the study, our results showed a large effect size (see appendix for Cohen values) regarding the increase of perceived pain intensity and believability in the presence of subtle trunk movement. Idle motion, in particular the movement reproducing human equilibrium, may be critical when studying emotion perception, as it favors the avatar's anthropomorphism and makes sure that it really conveys emotion.

Effect of trunk oscillation as pain behavior on participants' perception of the avatar's pain

In Experiment 1, significant increase of pain intensity scores from conditions 1° to 3° indicated that the wider the trunk movement was, the more intense the avatar's pain was judged, even in the condition with a voluntary exaggerated oscillation (3°). Pain believability rating was also widely improved when the avatar's trunk oscillated from 1° to 3° compared to the static (0°) and idle (0.67°) conditions. These results support the idea that when systematically detectable, the trunk rocking in this context may evoke a pain behavior.

Although this interpretation is yet to be verified, it was proposed that the trunk oscillation might reflect the bodily expression of a sustained pain. Overall, this work is concordant with seminal studies showing that the facial expression of a CG character may convey more emotion when combined with body cues [START_REF] Buisine | The Role of Body Postures in the Recognition of Emotions in Contextually Rich Scenarios[END_REF][START_REF] Clavel | Combining Facial and Postural Expressions of Emotions in a Virtual Character[END_REF][START_REF] De Melo | Influence of autonomic signals on perception of emotions in embodied agents[END_REF]. Moreover, contrary to pain intensity, pain believability did not significantly vary between conditions 1° to 3°, suggesting that it reached its peak at 1° of trunk oscillation.

Although partially different than those expected a priori, these results confirmed that the relation between pain believability and perceived pain intensity is not linear. It is not enough to make the CG character express more and more pain to make his pain seem more realistic and convincing.

Importantly, the present study goes even further by revealing that this assumed pain behavior may be easily manipulated to intensify the avatar's facial pain expression, as well as enhancing its accuracy. Indeed, in the static condition, the level of perceived pain in the avatar was much lower (mean score = 25.21, 95% CI [18.06,32.36]) than the level of the modeled pain expression (60 % of AU's maxima), while there was a relatively good accuracy when the pain expression was combined with a large trunk movement (e.g. 1° = 36.93, 95% CI [30.39, 43.47] ; 3° = mean score 55. 38,95% CI [48.02,62.74]). A similar underestimation of the avatar's pain intensity has also been observed in a recent study using the EEVEE platform [START_REF] Tessier | Toward dynamic pain expressions in avatars : Perceived realism and pain level of different action unit orders[END_REF]. Participants were presented with video clips displaying a sequential or synchronized dynamic facial pain expression of four different avatars modeled at the maxima of pain AUs (100 %). The mean score of pain intensity ranged from 47 to 57 (values from a 0-100 scale) according to the condition, which was half as less as the level of the modeled pain expression. Herein, the presence of trunk oscillation would enrich the avatar's pain expression, resulting in a better congruence between the subjective perception of its pain and the 'presumed' pain level of its facial expression. It would be interesting to explore whether other pain behaviors such as freezing [START_REF] Defrin | The evaluation of acute pain in individuals with cognitive impairment : A differential effect of the level of impairment[END_REF][START_REF] Feldt | The checklist of nonverbal pain indicators (CNPI)[END_REF] would induce similar results. If so, this would offer a varied range of pain behaviors while keeping a highcontrolled experimental set-up. It would also be interesting to investigate the influence of bodily expression combined with a neutral facial expression to determine whether people can infer pain from bodily expression alone or if facial expression is necessary. [START_REF] Atkinson | Emotion Perception from Dynamic and Static Body Expressions in Point-Light and Full-Light Displays[END_REF], showed that participants were able to detect emotion from points of light representing the shape of actors expressing emotions. We could then suppose that participants would be able to infer pain from body language alone, but this remains to be tested by adding a condition with a neutral expression.

Results for pupil size also led to an unexpected pattern, revealing a significant higher pupil dilation for the conditions 1° and 3°, and no change for the conditions 1.67° and 2°

(Figure 5). We proposed that the greater pupil dilation for the condition 3° would rather be an artifact resulting from the 'gap' created between this condition and the other pain behavior conditions. Indeed, the oscillation of 3° would contrast with the other conditions, whose the difference in amplitude was approximately the same. Condition 3° would have been perceived

implicitly as an 'isolated' event, provoking a greater pupil dilation each time this video was displayed. In other terms, pupil change would not reflect a greater emotional arousal due to the perception of pain in the avatar, but would be rather a response to its exaggerated trunk movement. This alternative explanation is in agreement with a study showing that pupil dilation may reflect not only emotional arousal, but also mechanisms underlying unconscious processing of surprising events (Alamia, Van Rullen, Pasqualotto, Mouraux & Zénon, 2019).

Consequently, it is assumed that the pupil responded more to the avatar's pain expression only when its trunk oscillated at 1° compared with the static condition. The higher pupil response may reflect a higher emotional arousal due to the enhanced realism consistency of the avatar's pain expression, i.e. between facial and bodily expression of pain. Although speculative, this may support the fact that this trunk amplitude is the optimal rocking level to combine with a facial pain expression modeled from AUs contracted at an intensity of 60 %.

Conclusion

This study complements the existing human-machine interaction literature aiming to disentangle how realism consistency between appearance and behavior of a CG character may affect the way we interact with it. Our results showed, for the first time, that both noncommunicative and emotional body motion combined with facial expression made the avatar's pain perception seem more accurate, and in some cases provoked a greater emotional arousal than facial pain expression with a still body. This is in line with the bulk of evidence demonstrating the potential of body language, when congruent with facial expression, to make avatars' emotions seem more realistic, and thus 'acceptable' enough to study underlying processes in humans [START_REF] Bogdanovych | What makes virtual agents believable?[END_REF]. Another novelty of the present work was to focus on a pain-expressing CG character. Pain expression is rarely explored in studies questioning avatar realism, which is paradoxical considering the growing number of publications in the medical and pain research domains and the added value of virtual reality for medical training, or the understanding and remediation of pain communication [START_REF] Brunet-Gouet | Editorial : Advances Computing for the Understanding and Remediation of Social Cognitive Disorders[END_REF][START_REF] Meugnot | The Contribution of New Technological Breakthroughs to the Neuroscientific Research of Pain Communication[END_REF]. Although several limits have been pointed out, this work may encourage researchers to pursue the endeavor of developing fully expressing CG character platform dedicated to cognitive science and health research.

Appendix -Post-hoc comparison tables

Figure 1 -

 1 Figure 1 -Left figure: Avatar at the onset of the clip (neutral expression). Right figure: Avatar expressing pain (AUs recruited at 60% of their maximal contraction)

Figure 2 :

 2 Figure 2: illustration of the trunk movement implemented on the avatar

  8). The data were epoched relative to the onset and end of the video clips (70 epochs x 2 tasks = 140 epochs of 4 seconds for each participant). The first second of each epoch was removed from the analysis to exclude the change related to the initial pupillary light response[START_REF] Mitz | Using pupil size and heart rate to infer affective states during behavioral neurophysiology and neuropsychology experiments[END_REF]. The average pupil size for each condition of trunk oscillation was then calculated between the 1-and 4-second marks. A change rate (%) in pupil size was calculated for the six conditions of trunk movement relative to the static condition (0°) using the following formula: (mean pupil size for given condition * 100 / mean pupil size for condition 0°) -100. A positive change indicated a larger pupil size in the condition of trunk movement than in the condition 0°. A negative change indicated a smaller pupil size in the trunk movement than in the condition 0°.All statistical analyses were conducted using SPSS. Analysis of variance (ANOVAs) with repeated-measures were performed on both Pain Intensity (0-100) and Pain Believability (%) with the condition of trunk oscillation (0°, 0.33°; 0.667°; 1°; 1.667; 2°; 3°) as withinsubject factor. A Greenhouse-Geisser correction was used when the sphericity assumption was violated (p < .05 for the Mauchly's test of sphericity). Alpha was set at .05 and in case of significance, pairwise comparisons using Bonferroni correction were applied. T-tests relative to 0 were performed on Pupil Size Change for each condition of trunk oscillation (0.33°;

  Pain Intensity task. The repeated-measures ANOVA revealed a large and significant effect of trunk movement on the evaluation of Pain Intensity [F(1.5,35) = 43.54, p < .001, η² = .30, ηp² = .65]. As illustrated in Figure 3, pairwise comparisons with Bonferroni correction indicated that the judgment of Pain Intensity increased significantly with the amplitude of trunk movement, except between the lower amplitude of oscillation (0.33°) compared to the condition without oscillation (0°) (p = 1.0).

Figure 3 -

 3 Figure 3 -Perceived Pain Intensity of facial pain expression according to the amplitude of trunk movement.Error bars correspond to the confidence intervals at 95% (using t table and corrected standard deviation because N < 30). All significance thresholds were adjusted using Bonferroni correction.

Figure 4 -

 4 Figure 4 -Pain Believability rating according to the amplitude of the oscillations Error bars correspond to the confidence intervals at 95% (using t table and corrected standard deviation because N < 30). All significance thresholds were adjusted using Bonferroni correction.

Figure 5 -

 5 Figure 5 -Mean of pupil size change, relative to 0 (%) according to the amplitude of trunk movement. Error bars correspond to the confidence intervals at 95% (using t table and corrected standard deviation because N < 30)

  South-Paris University (11 women, range = 18-27 years, mean age = 20.3 ± 2.65 years) through advertisements sent via South-Paris University's student mailing lists or a course program which proposed, on a voluntary basis, to obtain one extra mark on a test when participating to an experiment conducted in the laboratory. All participants had (or were corrected to) normal vision. This study was approved by the Ethics Committee of Paris-Saclay University (#CER-Paris-Saclay-2018-16) and participants gave written informed consent to take part in this study.

Figure 6 -

 6 Figure 6 -Proportion of detection according to the amplitude of trunk movement

Declarations of interest

None.


Pierre Tuffery

Sjoerd De Vries

The search of sequence variants using a constrained protein evolution simulation approach

Protein engineering or candidate therapeutic peptide optimization are processes in which the identification of relevant sequence variants is critical. Starting from one amino-acid sequence, the choice of the substitutions must meet the objective of not disrupting the structure of the protein, not impacting the main functional properties of the starting entity, while also meeting the condition to enhance some expected property such as thermal stability, resistance to degradation, ... Here, we introduce a new approach of sequence evolution that focuses on the objective of not disrupting the structure of the initial protein by embedding a point to point control on the preservation of the local structure at each position in the sequence. For 6 mini-proteins, we find that, starting from a single sequence, our simple approach intrinsically contains information about site-specific rate heterogeneity of substitution, and that it is able to reproduce sequence diversity as can be observed in the sequences available in the Uniref repository. We show that our approach is able to provide information about positions not to substitute and about substitutions not to perform at a given position to maintain structure integrity. Overall, our results demonstrate that point to point preservation of the local structure along a sequence is an important determinant of sequence evolution.

Introduction

Protein or peptide engineering or design, as well as therapeutic candidate peptide identification and optimization are processes in which the identification of relevant sequence variants is critical. Indeed, even not considering the insertion of un-natural amino-acids, the choice of substitutions must face the objective of not disrupting the structure of the protein, while meeting the condition to enhance some expected property such as stability, resistance to degradation, reduced immune response, among others. Ways to identify such substitutions include experimental techniques such as phage display [START_REF] Hoess | Protein design and phage display[END_REF], directed mutagenesis [START_REF] Bornscheuer | Improved biocatalysts by directed evolution and rational protein design[END_REF], and more recently deep mutation scanning [START_REF] Fowler | Deep mutational scanning: a new style of protein science[END_REF]. Due to progress in understanding the determinants of protein folding and increasing amount of data about related sequences, computational approaches to assist the design have also been developed. Among these, the simulation of sequence evolution along a phylogenetic tree has been the subject of intense efforts during the past decades [START_REF] Fletcher | INDELible: a flexible simulator of biological sequence evolution[END_REF][START_REF] Sipos | PhyloSim-Monte Carlo simulation of sequence evolution in the R statistical computing environment[END_REF][START_REF] Spielman | Pyvolve: a flexible Python module for simulating sequences along phylogenies[END_REF][START_REF] Low | COMPASS: the COMPletely Arbitrary Sequence Simulator[END_REF]. Such simulation techniques are now able to deal with different rates of substitution along sequences [START_REF] Shapiro | Choosing appropriate substitution models for the phylogenetic analysis of protein-coding sequences[END_REF], consider co-evolution for some [START_REF] Camenares | Simulating protein and nucleic acid sequence co-evolution[END_REF], and models to manage indels have emerged [START_REF] Fletcher | INDELible: a flexible simulator of biological sequence evolution[END_REF]. Apart from the fundamental goal of deciphering the rules underlying protein family sequence evolution, these approaches also have implications for enhanced sequence alignment for instance [START_REF] Anderson | SuiteMSA: visual tools for multiple sequence alignment comparison and molecular sequence simulation[END_REF][START_REF] Löytynoja | Accurate extension of multiple sequence alignments using a phylogeny-aware graph algorithm[END_REF][START_REF] Katoh | A simple method to control over-alignment in the MAFFT multiple sequence alignment program[END_REF]. These techniques based on the analysis of sequences alone can however reveal to be limited for protein engineering. In the context of sequence optimization, the search for sequence variants might escape the rules of natural evolution, and phylogenetic inference is not always possible in the context of de novo protein design, e.g. for the design of un-natural proteins. Finally, structural constraints, although implicitly contained along sequence evolution might require more direct consideration. Explicit account for structure has led to the field of computational protein design, in which the impact of candidate substitutions on the free energy of folding of the protein is usually questioned to drive the process [START_REF] Xia | Simulating protein evolution in sequence and structure space[END_REF]. Indeed, various techniques to question the relationship between the sequence space and the structure space, have been developed in the perspective of protein design or engineering. These include sophisticated protocols (see [START_REF] Gainza | Algorithms for protein design[END_REF], such as Rosetta-design [START_REF] Liu | RosettaDesign server for protein design[END_REF], protocols based on the dead end elimination theorem (e.g. [START_REF] Allouche | Computational protein design as an optimization problem[END_REF], protocols focusing accurately on the relationship between evolution and side chain-packing [START_REF] Perron | Modeling Structural Constraints on Protein Evolution via Side-Chain Conformational States[END_REF], coarse grained models (e.g. [START_REF] Grahnen | Biophysical and structural considerations for protein sequence evolution[END_REF], to cite some, or a combination of the simulation of protein sequence evolution under the constraint of explicit structural context [START_REF] Grahnen | CASS: Protein sequence simulation with explicit genotype-phenotype mapping[END_REF]. Of note, 3D based-techniques have led to successful applications in protein design [START_REF] Bhardwaj | Accurate de novo design of hyperstable constrained peptides[END_REF][START_REF] Huang | The coming of age of de novo protein design[END_REF][START_REF] Dou | De novo design of a fluorescence-activating β-barrel[END_REF].

In the context of sequence optimization of therapeutic peptide/protein candidates, a primary objective, beyond the goals of enhancing peptide stability and bioactivity, is the preservation of the 3D conformation among the sequence variants, not to disrupt the geometry of the functional interactions with a target. Thus means to identify which substitutions are likely to disrupt or to preserve those interactions are highly desirable. In particular, it is well-known that some sites of the protein sequence are much more sensitive for this than others. For example, in protein folding, a few key residues force the chain to adopt a rudimentary native-like architecture [START_REF] Dobson | Protein folding and misfolding[END_REF]. In sequence evolution, the phylogenetic tree is characterized by a few tree-invariant positions -detectable by Evolutionary Tracing [START_REF] Lichtarge | An evolutionary trace method defines binding surfaces common to protein families[END_REF], and tree determinants (25) that correspond to functionally important residues in the structure. Therefore, it is essential to identify those key residues and delineate the sequence space that preserves the 3D conformation in their vicinity. Although protein design techniques coupled to molecular simulations techniques can fulfill such a goal, such protocols are hardly tractable to test the complete space of sequences that rapidly growths with peptide size.

Here, we introduce a mid term strategy which does not rely on sophisticated phylogenetic inference nor on explicit structure consideration. We use the concept of structural alphabet [START_REF] Camproux | A hidden markov model derived structural alphabet for proteins[END_REF], a discrete model of local structure, to constrain a very simple protocol for sequence evolution. Our previous works have demonstrated the identification and use of the concept of structural alphabet, and the fact that divergent profiles tend to correspond to divergent structures, as illustrated by a fragment search strategy exploiting this trend [START_REF] Shen | Detecting protein candidate fragments using a structural alphabet profile comparison approach[END_REF]. This is the first time that our structural alphabet is applied to the effect of mutations of the same structure. We show that it provides, at a limited computational cost, sequence sampling consistent to that observed in the subset of the natural sequences identified and that it provides information about critical positions, ie. conserved positions not to substitute or position specific disruptive amino-acids substitutions.

Theory

Structural alphabet (SA)

A structural alphabet (SA) can be considered as a generalized secondary structure in which each region of a protein can be associated to one specific conformation, or "SA-letter" of the alphabet. Thus, a structure can be transposed into an optimal string of SA-letters expected to describes accurately enough the series of the local conformations -two different strings should be able to distinguish between small conformational changes and hopefuly provide a basis to regenerate the 3D conformation. Here, we use a structural alphabet made of 27 SA-letters that have been identified using Hidden Markov Models [START_REF] Camproux | A hidden markov model derived structural alphabet for proteins[END_REF]. Each SA-letter corresponds to a fragment of four amino-acids of the protein, and is described by four distance descriptors and the SA-letters overlap by three amino-acids. Thus, the conformation of a protein of size L amino-acids can be described by a string a SA-letters of size L-3.

Given the amino-acid sequence of L amino-acids, it is possible to predict the probability that each series of 4 amino acid is associated with each of the 27 SA-letters, i.e. a SA-profile of size 27(L-3) values. In such profile, each of the 27 probabilities associated to position l correspond to estimates of p(SA/aa), i.e. probabilities of the SA letters given information of the amino-acids. Briefly, for a given protein sequence, a sequence profile is first generated using PSI-BLAST. It is then used to predict an SA-profile. This prediction is performed by a Support Vector Machine (SVM) that was trained to reproduce the SA-profiles extracted from known structures. Here, we use a 2012' update of the SVM predictor that was described in [START_REF] Maupetit | A fast method for large-scale De Novo peptide and miniprotein structure prediction[END_REF], learnt over a collection of proteins of size more than 80 amino acids. Of note, those profiles have been demonstrated effective for the de novo modeling of peptides, using different strategies to exploit the information they contain [START_REF] Lamiable | PEP-FOLD3: faster de novo structure prediction for linear peptides in solution and in complex[END_REF].

Interestingly, the SA-profiles make a link between sequence and structure, and their variations are likely to reflect how the substitutions in the amino-acid sequence can impact the 3D structure of the protein. They contain information specific of the structure, the profiles of two proteins of dissimilar folds being dissimilar and the profiles of tow proteins of similar fold being similar. We have shown in a previous study [START_REF] Shen | Detecting protein candidate fragments using a structural alphabet profile comparison approach[END_REF] that using the Jensen Shannon divergence (JS) as a measure of profile dissimilarity, it is possible to identify, for low values of JS, protein fragments that have similar conformations. Our purpose here is to use SA-profile divergence as a control over the preservation of the structure of the target after an amino-acid substitution, to ensure a position by position preservation of the local conformation associated to the initial amino-acid sequence. For the purpose of illustration, Figure 1 shows the relationship between the JS and the RMSD for a fragment of 4 amino acids taken from the structure of the GCC-box binding domain (PDB(30): 3GCC), positions 151-155 aligned with random fragments of a collection of proteins. One clearly sees that small JS value (resp. large JS values) tend to be associated with small RMSD values (resp. large RMSD values). Note that there are a few cases for which the RMSD is large while the JS is low, which is in agreement with the well known observation that small fragment with identical sequences can adopt different conformations [START_REF] Sternberg | Local protein sequence similarity does not imply a structural relationship[END_REF]. Conversely some fragments come with a large JS value (different sequences) and small RMSD values, again consistent with previous observations. Importantly, in the present study, we do not consider unrelated fragments, but we make evolve one specific sequence, which is expected to prevent the occurrence of such singular behaviors. Figure 2 illustrates how a SA-profile can be affected by one single mutation. Here, we present the impact of the A53T mutation of the α-Synuclein, reported to impact the local conformation around position 53 [START_REF] Kang | The A53T mutation is key in defining the differences in the aggregation kinetics of human and mouse α-synuclein[END_REF]. One clearly sees that just one single substitution results in large differences in the SA-profile around position 53, not affecting, as could be expected other parts of the profile (JS values of more than 0.13). Note that in our SA model, the probability distributions correspond in fact to overlapping fragments of 4 amino acids. Thus, one amino acid contributes to 4 consecutive probability distributions.

Methods

Sequence evolution

Figure 3 presents an overview of the procedure. The protocol to simulate sequence evolution is similar to a Monte Carlo procedure. Starting from the sequence of a target protein, it will generate substitution events that are accepted or not, simulating a random walk in the sequence space. For each substituted sequence, a SA-profile is predicted and compared to the SA-profile generated with the initial sequence. Acceptation of substitutions is performed according to the expected structural divergence as measured by JS.

The core simulation of sequence evolution is made according to the scheme proposed by PSeqGen [START_REF] Grassly | PSeq-Gen: an application for the Monte Carlo simulation of protein sequence evolution along phylogenetic trees[END_REF] or CS-PseqGen [START_REF] Tufféry | CS-PSeq-Gen: simulating the evolution of protein sequence under constraints[END_REF]. The JTT substitution matrix is used to allow for different rates of substitutions depending on the type of amino acids. Since we are interested in the substitution events disregarding any molecular clock or any phylogenetic model, evolutionary time is considered as infinite, meaning each step results in selecting a substitution, similarly to a Monte Carlo inspired process. Site-specific rate heterogeneity is not considered either, since in the context of protein design, it is largely unknown. Substitution acceptance is controlled as detailed below. In total, simulations are defined by an initial sequence and a number of substitutions to generate.

Acceptance or reject of the substitutions is controlled by the estimation of the expected impact of a substitution on the preservation of the local structure of the target, i.e. on the dissimilarity of the SAprofiles before and after substitution, using JS as follows:

JS(SAPref, SAPcur) < JScut-off: accept, reject otherwise where SAPref corresponds to the SA-profile of the initial sequence, SAPcur to the SA-profile of the current sequence and the JScut-off value corresponds to the value above which the divergence is considered as too large to guarantee structure preservation. Varying the value of this cut-off can be assimilated to varying the temperature of the Monte-carlo. SA-profiles associated with an amino-acid sequence of size L consist in L-3 probability distributions, each of 27 values. For each pair of distributions belonging to SAPref l , SAPcur l , , where l varies from 1 to L-3, we use JS(l) as a mesure of the dissimilarity. As mentioned above, a particular position l in the sequence contributes to the SA-profiles of four positions (l .. l+3), which requires to check the JS at those four positions. For sake of simplicity, we measure the JS for all of the L-3 positions, and to combine the L-3 JS values, we consider:

JS(SAPref, SAPcur) = max(JS(l), l = 1 .. L-3)
where the maximum stands to ensure that, given a cutoff value, for no profile the deviation is more than this value.

Of note, the observed variability over independent simulations is, in our experience, rather weak and is not discussed any further here. Also note that for each substitution event, we perform a full SVM prediction which requires itself a psi-blast and thus simulations of 3000 substitution events take between one and two days each on a standard workstation.

Amino-acid profile divergence

Once simulations performed, one needs means to analyze the diversity of the sequence generated. To estimate the number of amino acids occurring at a given position l in the sequence, we use the number of equivalent amino-acids defined as:

Neq(l) = 1/pmax(l)
where pmax(l) corresponds to the maximum of the occurrence probabilities associated with the 20 amino-acids at position l in the sequence. These values are averaged over the L positions of the sequence to get an estimate of the diversity over all positions.

The comparison of the two profiles made of L distributions of 20 probabilities is a difficult issue. Here, we simply consider the correlation coefficient of the two vectors made of the L . 20 values, which we found more intuitive than criteria such as the average JS or the average dot product over the L positions.

Test sets

To assess the performance of our procedure, we have chosen 6 peptides and mini-proteins of known structure, not intrinsically disordered, and of varied topologies. Table 1 provides details on their size, secondary structure, and sequence. Such small sized proteins have the advantage of not belonging to the range of size used to learn the SVM. We emphasize, however, that although we assess the protocol on such small proteins, the procedure is very general and applicable to larger proteins, as illustrated Figure 2 for the α-Synuclein, a protein of 140 amino acids.

As a reference for "naturally observed" sequence variation, we have considered the sequence variants made of the homologs of these targets found either in the Protein Data Bank, or in the Uniprot repository. For the later case, since our aim here is to estimate the acceptable substitutions in terms of amino-acids, we have build profiles considering entries of the Uniref 90 subset, to discard too identical sequences. The procedure was as follows:

1/ search uniref90 using blastp (35) 2/ retrieve the entries of the hits 3/ perform a multiple alignment of the hits using clustal omega (36) 4/ identify the region of the alignment matching the query and extract the corresponding profiles Name PDB s2 L Seq.

Results and discussion

Sequence divergence depends on the JS cut-off value

Figure 4 depicts, for the transactivation domain of CRE-BP1/ATF-2 (PDB: 1bhi), the evolution of the JS and of the sequence identity as a function of the number of substitution events. As could be expected, one observes a fluctuation of JS between low and large values (top inset, red) whereas JS of the accepted sequences remains under the limit imposed of 0.12 (top inset, black). Interestingly, one also observes that the sequence identity of the sequences generated, when compared to the initial sequence, rapidly decreases down to values close to only 10-20% (black), which corresponds to the twilight zone in terms of homology. Small discrepancies are observed between the sequences accepted and those generated, on average, as can be observed for instance around step 200. A control simulation using the same procedure, but not applying any JS constraint shows that rapidly, sequence diverge in a random manner (gray).

Varying the value of the JS cut-off impacts, as expected, the degree of divergence of the sequences. Indeed, lower value of JS cut-off results (JS cut-off 0.08, green) in exploring sequences with higher sequence identities, given that a more strict control on the impact of substitutions on the predicted local conformations results in rejecting more substitutions.

Similar behaviors are observed for the other targets of the test set, and are summarized in Table 1 that reports, for all targets, the lowest sequence identity reached during the simulations. Overall, the results show that the procedure is able to reach explore sequences diverging down to close to 10-20 % sequence identity. Figure 5B-D present logos corresponding to subsets of sequences obtained from a simulation of 3000 events driven by a JS cut-off of 0.12, using as second step filters JS cut-off values of 0.12, 0.10 and 0.08, respectively. The size of the subsets is of 892, 256 and 29 sequences, respectively. As a reference, Figure 5A presents the results of an unconstrained simulation of 892 events, i.e. a size identical to the subset depicted Figure 5B. Figures 5E and 5F correspond to the logos obtained from the PDB (13 sequences only) and uniref90 (323 sequences), using blastp (see methods).

Constrained simulated sequence evolution intrinsically embeds site-specific rate heterogeneity of substitution

A first observation is that compared to the unconstrained simulation , that leads to a rather flat profile, the logo profiles differ largely. In Figure 5B, positions 6, 14, 18 and 31, associated with amino acids P, C, F and H, respectively, appear more conserved than other positions (Neq < 1.3). Interestingly, considering lower JS cut-off values (Figure 5C-D), the number of positions with Neq < 1.3 increases to encompass progressively also positions 7, 19, 22. According to the design of our procedure, this suggests that these amino acids are probably critical to maintain the local conformation at these positions. Looking at the structure of 1bhi (Figure 6), one notes that these positions are located at the extremities of secondary structure elements. Positions 9, 14, 27 and 31 correspond to the cysteines and histidines involved in the coordination of the zinc ion (not present in the structure).

Looking at the profile generated using homologues of the PDB (Figure 5E), only positions 9, 14, 18, 20, 22, 27 and 31 are associated with Neq values less than 1.3. Despite the weak number of sequences, those positions match rather well the results of the simulation. Note that, since we could not identify homologues with known structure for all 6 cases, we do not discuss the PDB profiles for the other targets.

Looking at the distributions obtained for uniref90, one observes that not only positions 6-7, 14, 18, 22 and 31 are conserved, but in fact positions 6-7, 9, 11-14, 16-18, 21-24, 27, and 29-37. Apart from C-ter residues 32-37 that are involved in the functional interaction with kinases, as reported by ELM [START_REF] Dinkel | ELM-the database of eukaryotic linear motifs[END_REF], those residues correspond to amino acids located at the interface between the secondary structure elements, most being involved in long range interactions stabilizing the overall conformations. Hence, the observed heterogeneity obtained from simulations constrained by JS seems to effectively contain information about key residues that are required for the preservation of the secondary structure elements.

Sequence divergence is target dependent.

We now consider all proteins of the test set. Table 2 reports for all targets the number of sequences accepted for different JS cut-off values, and the average Neq values of the profiles. Note that all simulations have been performed using JS of 0.12. For JS values below 0.12, the numbers correspond to the subsets of the simulated sequences matching the JS condition. It is striking that the numbers depend on each target, meaning target sequence/topology has an impact on the simulated evolution. Table2: #: Number of sequences accepted at the given JS value, and the associated Neq. Simulations of 3000 substitution events were performed using JS=0.12. Uniref90: corresponding values observed among homolog sequences of uniref90. JSu: JS values for which simulated Neq is identical to that of uniref90. JS2: JS cut-off values used for the second stage simulations (section 4.4). JSu2: JSu for the second stage simulations. r: Correlation coefficient between second stage and uniref90 profiles.

If, for a JS-cutoff value of 0.12, the observed variation in the number of sequences is rather limited (acceptance rate of substitution events between 22 and 43 %), larger variations are observed for lower JS cut-off values. For instance for 0.08, the number of sequence varies between 405 and 27 (acceptance rate between 13 and 0.9%), and the Neq values vary between 2.87 and 1.75, without any obvious relation to the number of accepted substitutions.

Finally, if we consider the values of Neq obtained for the uniprot90 homologs, they look rather small (between 1.05 and 2.19) compared to the values obtained from the simulations. For the simulations, values of Neq similar to those observed for uniprot90 are obtained for rather low values of JS -between 0.03 and 0.12 (JSu), but again one observes a strong dependence on the specifics of each particular case.

JS controlled sequence evolution can mimmick sequence divergence observed in natural sequences.

In order to get a better agreement between the simulations and the observations of uniref90, we have performed second stage simulations using JS values for the Neq were close to that of uniref90 (Table 2, JSu2 values). Note that these new cut-off values guide the search into new regions of the sequence space, and thus, we do not expect the Neq values to be preserved, thus JS2 values were chosen slightly larger than the JSu values. A first observation about these simulations is that indeed, the acceptance rate is, as expected much larger for lower JS values (between 18 and 37% -see Table2) than that observed for the initial simulations, which confirms a better sampling of the sequence space for the targeted divergence, the minimum number of accepted substitution events being of 538, when it was below 100, for the first series of simulations. The minimal sequence identities for the accepted sequence were of 23, 37, 32, 20, 31 and 25% for 1bhi, 1by0, 1dv0 2k9d, 3e21 and 1wr3, respectively.

Figure 7 shows for each of the 6 targets, the logos obtained for the JSu2 values, i.e. for average Neq values similar to those obtained for uniref90. There is a good visual agreement between the logos of the simulations and uniref90, and the r values vary between 0.63 and 0.77 (see Table 2), which is highly significant (p < 10 -10 ). Note however that the number of sequences, even if much larger than that of the first series of simulations, remains low, except for 2k9d and 3e21 for which the r values, even if a bit lower remain very significant. Of note, this low number of accepted sequences is also conditioned by the weak Neq values observed for uniref90. We have further verified that for series of 1000 tests permuting the series of 20 probabilities for one of the profiles, the distribution of r values is close to 0-0.1 on average. Overall, this suggests that our simple procedure is able to mimmick accurately enough the sequence fluctuations as observed in uniref90.

Looking at the positions highly conserved (Neq < 1.1) it is striking that most of them are located at the termination of secondary structures, which suggest that the approach is probably sensitive enough to detect that these regions are critical for structure preservation. This is even more true when looking at residues conserved for slightly larger values of JS cut-offs (magenta in Figure 7). Overall, the positions highly conserved identified by our protocol correspond to 31 over 55 (56 %) of the positions highly conserved observed in the uniref90 distributions. As illustrated upper for the 1bhi case, constraints on sequence can occur due to several reasons not related to the preservation of the local structure.

Finally, another outcome of the approach is that it also provides for each position in the sequence, information about the substitutions that have been rejected. Figure 7 also presents, for all 6 cases the logos of the amino acids that were not accepted by the JS cut-off, and were excluded at least once with a JS value of more than 0.3 (chosen large enough compared to the JS cut-off values to get confident enough a local conformational change is expected). It is interesting to compare these amino acids to those present in the uniref sequences. For instance, at 1by0 conserved position 22, D was rejected several times by our protocol, whereas E was preserved. D is not present either in the uniref profile. At position 10 of 1dv0, the uniref profile has a conserved position with R, when our simulations report occurrences of R,K and M, and A, L, T and V where rejected. For 3e21 position 10, the uniref contains occurrences of D, N, S, E, accepted amino acids of the simulations contain D, N, E and G, A, H, R, S, T, Y were rejected. Another such similar example is observed at 1wr3 position 28. In few cases, however, some rejected amino acids are observed in the uniref profiles. Such a situation occurs for 2k9d position 10, where the simulations accepted R, H, K, uniref contains occurrences of N, H, K, R and rejected amino-acids include N. Several causes can explain such events, including epistatis effects (occurrence of N conditioned by the occurrence of another amino-acid close in the sequence to preserve local structure), or insufficient accuracy of the SVM prediction for some sequence motifs, to mention two of them. Such cases are however rare, and concern mostly amino acids occurring with weak frequencies in the uniref90 profiles. Over the 6 targets of our test set, the average cumulated probabilities of such amino acids if of only 6.7%.

Discussion and perspectives

In the present study, we have introduced a new paradigm to simulate protein sequence evolution. It builds on the concept of structural alphabet -a kind a generalized secondary structure -used as a means to constrain a very basic procedure to sequence evolution. The underlying hypothesis is that the limitation of the divergence between two structural alphabet profiles ensures the point to point preservation of the local conformation along the sequence. Surprizingly, we find that this very simple procedure is able, for the mini-proteins of our test set, to result in generating sequences which significantly fit sequence diversity as observed among uniref90. Indeed, we find that our procedure is able to grab the specifics of the heterogeneity of rate among sites, and of the nature of the amino-acids at those sites. This being only induced by the constraint on the divergence between the structural alphabet profiles implies that the point to point control of the preservation of the local structure of a target is an important determinant of sequence evolution. To the best of our knowledge, our procedure is the first procedure to directly address this question, and this clearly opens the door to further investigations on the amount by which sequence evolution is controlled by such local constraints versus longer range interactions. One has however to consider carefully the present results. For one part, we have tested so far a limited number of targets, and in addition, some optimizations can already be foreseen. A special point is about the heterogeneity in the JS cut-off values depending on the targets. Probably, it could be of interest to study the convergence of our simulations for a larger number of cycles of simulations -here, only two were performed. A comparison and/or combination with the results of simulation procedure along a phylogenetic tree could also reveal informative.

It is difficult to compare the results of our protocol to those of former approaches described to measure the impact of mutations, as their objectives differ largely. The aim of the approaches developed so far has usually been to predict the impact of substitutions on protein stability. This includes local conformation preservation considerations, but also longer range effects involved in the stabilization of the structure, such as the interactions between secondary structure elements, or more complex effects.

Here our main focus is only about the preservation of the local structure, disregarding longer range effect. It is thus expected that our results are more specific than those of other approaches. In addition, most of former approaches, such as Backrub [START_REF] Lauck | RosettaBackrub-a web server for flexible backbone protein structure modeling and design[END_REF], Rosetta-Design [START_REF] Liu | RosettaDesign server for protein design[END_REF] or POP-Music (39) require the 3D coordinates of the protein while we consider here only the sequence. We discuss briefly our results in the light of the results obtained by the INPS server (40) -based on sequence and the very recent SAAMBE-3D server (41) -based on structure. Considering for instance the 1by0 target, both INPS and SAAMBE-3D predict that the subsitution E22D is slightly destabilizing, while it is not accepted by our protocol, nor observed in unriref90, which seems consistent. At position 10 of 1dv0, SAAMBE-3D predicts all substitutions to be destabilizing -which is consistent the full preservation of K in uniref90, while INPS predicts that substitutions of R into A, T, K or M would decrease the stability of the protein and that substitutions into L or V would increase the stability. Our results suggest that substitutions of K into R or M are compatible with local structure preservation, occurrences of A, L, T or V being never observed. For 2K9D position 10, the same kind of observations can be done. R10H, R10N, R10K substitutions are all predicted to decrease the stability by both INPS and SAAMBE-3D, while our procedure suggests that R10K and R10H preserve local structure, and all are observed in uniref90. In fact, it is fully possible that some substitutions that destabilize the global structure do not affect the local structure and are still compatible with the preservation of the local structure. It is also to be noted that servers such as INPS or SAAMBE-3D consider the impact of single mutations, all other amino acids being preserved, while our procedure makes the complete sequence evolve, and thus, some substitutions can be accepted conditional to previous events in the neighborhood of a site. We also recall that our procedure focusing on the preservation of 3D-structure comes with no quantification of the impact of the substitutions on stability or on function.

Indeed, a motivation for developing this procedure was the perspective to assist sequence optimization, as can be needed for the development of a candidate therapeutic peptide or mini-protein. Two questions can be posed in such context. The first is to identify which positions in the sequence should not be modified to ensure the preservation of the conformation of the candidate. To this respect, it is interesting that the observed heterogeneity of substitution rates over positions obtained from our simulations make it possible, from a single sequence, to identity positions at which substitutions seem risky. A second is to propose, for positions at which we find some diversity, which substitutions are likely. Numerous protocols and prediction approaches have been developed to this aim (see for instance [START_REF] Steinbrecher | Predicting the effect of amino acid single-point mutations on protein stability-large-scale validation of MD-based relative free energy calculations[END_REF][START_REF] Quan | STRUM: structure-based prediction of protein stability changes upon single-point mutation[END_REF]. Here, a particularity of our approach, which is not quantitative, is probably to return information about which residue substitutions were not accepted at a given position. Although our results clearly show that rejected substitutions depend on the context of the neighboring residues since false rejections were observed at a frequency of close to 7%, we put emphasis on the fact that our procedure is based on single mutation events only, whereas more sophisticated ways to simulate evolution could be setup, for instance considering co-evolution. All together, we however hope that our procedure of sequence evolution under the constraint of the point to point preservation of local structure already meets the objective of assisting sequence optimization.

Funding: INSERM U1133 recurrent funding

Acknowledgments: The authors wish to thank J. Chomilier and Ph. Derreumaux for useful discussions. Right: Conserved (Neq < 1.1) residues as obtained from the simulations. Cyan residues depict residues that are still conserved in the subset of sequences using JS + 0.02, when possible (i.e. except for 3e21)

Figure 5

 5 Figure5presents, for 1bhi, a logo representation of the distribution of the 20 amino-acids per site. Figure5B-D present logos corresponding to subsets of sequences obtained from a simulation of 3000 events driven by a JS cut-off of 0.12, using as second step filters JS cut-off values of 0.12, 0.10 and 0.08, respectively. The size of the subsets is of 892, 256 and 29 sequences, respectively. As a reference, Figure5Apresents the results of an unconstrained simulation of 892 events, i.e. a size identical to the subset depicted Figure5B. Figures 5E and 5F correspond to the logos obtained from the PDB (13 sequences only) and uniref90 (323 sequences), using blastp (see methods).

Figure 1 :

 1 Figure 1: JS as a function of RMSD, using as seed the fragment 151-154 -4 amino acids -of the GCCbox binding domain (PDB: 3GCC) with all fragments of 77 small proteins.

Figure 2 :

 2 Figure 2: SA-profiles of the α-Synuclein, wild type (top) and A53T mutation (bottom). Each column corresponds to the probability distribution of the 27 letters. SA-letters are sorted from most helical (red) to most extended conformations (green). Details of the probability distributions around position 53 are provided as insets. Of note, SA letters correspond to fragments of 4 amino-acids. Thus the A53 position is associated with HGVA, GVAT VATV and ATVA fragments, facing the H, G, V and A columns, respectively.

Figure 3 :Figure 4 :Figure 5 :

 345 Figure 3: Flowchart of a simulation

Figure 6 :

 6 Figure 6: Example of 1bhi. Simulations using JS=0.12. The amino acids at the conserved positions (Neq < 1.3) are depicted in green. Amino acids conserved in the uniref90 profiles are depicted in blue.

Figure 7 :

 7 Figure 7

Table 1 :

 1 Name: the name of the protein. PDB: PDB identifier of the structure. s2: secondary structure topology. L: size (amino-acids). Seq.: amino acid sequence.

	Transactivation domain	1bhi	α2β 38 MSDDKPFLCTAPGCGQRFTNEDHLAVHKHKHEMTLKFG
	of CRE-BP1/ATF-2
	N-terminal leucine-	1by0	α	27 RKKLEELERDLRKLKKKIKKLEEDNPW
	repeat region of
	hepatitis delta antigen
	C-terminal UBA	1dv0	α3	47 GSQEKEAIERLKALGFPESLVIQAYFACEKNENLAANFLLSQNFDDE
	domain of the human
	homologue of
	RAD23A
	Bomain X of measle	2k9d	α3	44 VIRSIIKSSRLEEDRKRYLMTLLDDIKGANDLAKFHQMLVKIIM
	phosphoprotein
	FAF-1 UBA Domain	3e21	α3	45 GSMDREMILADFQACTGIENIDEAITLLEQNNWDLVAAINGVIPQ
	first WW domain of	1wr3	β3	36 GSPPLPPGWEEKVDNLGRTYYVNHNNRSTQWHRPSL
	Nedd4-2