Combining electrodermal activity analysis and dynamic causal modeling to investigate the visual-odor multimodal integration during face perception

The processing of facial expressions is a fundamental mechanism for perceiving others' intentions and emotions [1]. In real-life situations, this process is not limited to the visual system alone but is influenced by contextual information from other sensory channels [2]. Olfactory stimuli have been found to play a crucial role in modulating the hedonic perception of faces [1, 3]. Previous studies have shown that the emotional valence conveyed by contextual odors can affect the recognition of facial expressions and subjective ratings of faces [49].

These behavioral responses are accompanied by physiological changes, as evidenced by electroencephalographic (EEG) event-related potential (ERP) studies. These studies have reported an effect of the valence of the odors at both early sensory stages (P1/N1, N170, and vertex positivity potential (VPP)) [8, 1012] and later cognitive stages (late positive potential) of face processing [1315]. However, ERP waveforms reflect the overall activation of the face-perception system, an ensemble of interconnected regions which categorize visual stimuli as faces by analyzing various factors such as expression, gender, and identity [16]. The close association between olfactory and visual areas has been established in previous research [17, 18]. Therefore, investigating the effect of emotional odors on the interactions among the areas of the face-perception system could provide valuable insights into the integration mechanism between olfaction and vision. Such investigations can enhance our understanding of how sensory cues interact to shape our perception of others' facial expressions and emotional states.

The standard ERP analysis based on components' amplitude and latency can not provide a sufficient level of detail to investigate brain dynamics at the network level. Nevertheless, ERPs can be combined with more sophisticated techniques to provide a window on the cortical sources' dynamics underlying face processing. In this context, EEG connectivity analysis allows to investigate the interaction among neuronal assemblies [19]. Particularly, effective connectivity estimates the direct and directional (i.e. causal) influence that a source exerts over another [19]. Among the various effective connectivity approaches [19, 20], dynamic causal modeling (DCM) is a powerful model-based technique designed to test for the effects of experimental factors on the interactions among regions of a specified brain network, starting from observed electrophysiological or functional imaging responses [21, 22]. Several studies applied DCM to reveal the dynamics among sources of the face-perception network either through fMRI [2327], MEG [28, 29], or intracranial EEG [30] recordings. In this light, DCM can be combined with the information provided by visual ERPs to investigate the modulatory effect of hedonic contextual odors on the cortico-cortical interactions among face-processing areas.

A potential factor of interest in the analysis of the face-odor multimodal integration concerns the arousal elicited by faces. Specifically, within an event-related paradigm, several repetitions of the stimulus are presented over time, and it could be hypothesized that not all the faces are able to elicit the same affective response due to the subjective saliency of facial features [3133]. This may affect the signal-to-noise ratio (SNR) of the observed ERPs, and potentially reduce or even obscure the modulatory effect of contextual odors on the connectivity. However, although arousal has already been shown to influence face-evoked ERPs [34, 35], its effects when faces are presented with concomitant olfactory stimuli are still poorly investigated. Arousal has a crucial effect on the processing of external inputs. Motivationally relevant stimuli, such as the perception of a face with specific intrinsic features, are able to generate a transitory and automatic enhancement of arousal [36, 37], that entails a prioritized processing of the stimulus in the visual stream [35]. This mechanism is associated with a series of physiological changes, including a top-down (i.e. endogenous) affective influence on sensory gain control [35], an increase in the amplitude of frequency-specific brain oscillations [32, 38] and in the magnitude of ERP components [34, 35]. In this light, testing for the actual occurrence of enhanced arousal is crucial to ensure that a face has been successfully perceived, thus maximizing the SNR of ERP responses and the consequent effect size of odors' modulation.

Modeling the effects of arousal at both the ERP and DCM connectivity levels requires knowing in advance whether a stimulus has the property of eliciting it. However, such a prior knowledge is far from being easily identified, since particular features that could characterize a perceived stimulus as relevant (e.g. facial identity, eye gaze, expression) may vary on a subjective basis. Besides brain activity, arousal is known to influence autonomic nervous system (ANS) dynamics. Particularly, states of high arousal are associated with an increase of peripheral sympathetic activity [39]. In this light, electrodermal activity (EDA) can be exploited as an objective means to identify the occurrence of enhanced arousal to a given stimulus, possibly related to intrinsic facial features or to emotional expressions [4043]. EDA is comprised of a slow-varying tonic component overimposed to a fast-varying phasic component. Particularly, the latter is represented by a series of stimulus-evoked skin conductance responses (SCRs), whose elicitation is driven by peripheral sympathetic neural bursts: i.e. the sudomotor nerve activity (SMNA) [39]. Accordingly, sympathetic responses to faces observed from SMNA can be adopted as a reliable marker of enhanced arousal.

In this work, we investigate the effects of hedonic olfactory stimuli and arousal enhancement on face perception, as measured by ERP components and source effective connectivity. Specifically, we hypothesize that odors' valence modulates the strength of specific pathways within the face-perception network, and that arousal evoked by the saliency of faces could play a role on such a multimodal sensory integration. To this aim, we acquired the EDA and EEG signals from 22 healthy volunteers performing a passive visual stimulation task with neutral faces and background pleasant, neutral and unpleasant odors. We propose a novel methodological approach based on the convex-optimization-based EDA (cvxEDA) framework [44] to identify face-evoked peripheral sympathetic responses and characterize the stimuli according to their property of eliciting arousal. The outcome of this classification procedure is then adopted to model the effects of odors' valence and arousal as between-trial factors on the visual ERP components evoked by faces. We then exploit such ERPs to carry out a DCM and parametric empirical bayes (PEB) analysis of odors and arousal on the connectivity among cortical sources found to be activated by the task. Particularly, PEB allows to build a hierarchical analysis of effective connectivity, where single-subject estimates about neuronal parameters of interest (e.g. connections' strength) are treated as stochastic effects at the group level [45]. To validate our aforementioned hypothesis on the role of enhanced arousal to faces, we aim to assess: (i) whether such mechanism is effectively reflected by the observed ERPs, and (ii) the modulation operated by odors on cortico-cortical interactions at different levels of arousal.

2.1. Participants

Twenty-two healthy volunteers (age $27 \pm 3, $ females) were enrolled in the study. Volunteers did not report any history of neurological and cardiovascular diseases, anosmia or being tested positive to COVID-19 over the past 6 months. Volunteers were not allowed to have any food nor drink in the 30 min preceding the experiment. The study was conducted according with the guidelines of the Declaration of Helsinki and approved by the Bioethics Committee of the University of Pisa Review No. 14/2019, 3 May 2019. All participants gave written informed consent to participate in the study.

2.2. Olfactory stimuli

We selected three different odorants: i.e. banana (isoamyl acetate; $CH_3COOCH_2CH_2CH(CH_3)_2$), isovaleric acid ($(CH_3)_2CHCH_2COOH$), and n-butanol ($CH_3CH_2CH_2CH_2OH$). We chose such odorants to convey positive (banana), negative (isovaleric acid), and neutral (n-butanol) valence according to previous literature results [46, 47]. For each subject, we prepared 1 ml isointense solutions diluting pure odorant substances in distilled water according to the ratios of 1/20 (banana, n-butanol) and 1/10 (isovaleric acid) respectively. Odorants were delivered through an in-house built computer-controlled olfactometer at a flow rate of 60 ml min−1.

2.3. Olfactometer device

The 4-channels computer-controlled olfactometer used herein is composed by (i) a pressure regulator to set air pressure at 7 Bar, (ii) four stainless-steel containers (50 ml) equipped with o-rings, stainless steel caps and clamping rings to ensure a gas-tight closure, (iii) 10 low dead-volume three-way solenoid valves (Parker Hannifin, Italy), (iv) a digital flow meter (Honeywell, Italy), and (v) a disposable nasal cannula. An Arduino in-house code controlled the solenoid valves allowing them being opened and closed through a well-defined sequence of actions. The software allows the delivery of pure clean air (i.e. all valves are opened) or odors kept at room temperature within the containers. Components were connected to each other using polytetrafluoroethylene (PTFE) fittings and tubings (internal diameter of 0.3 mm) to reduce the olfactometer dead-volume up to 1 ml. The nasal cannula is connected to the olfactometer though a 3 m long PTFE line. The olfactometer showed a negligible memory effect and a low background emissions of chemicals in the air/odors mainstream as determined by thermal desorption coupled to gas-chromatography and mass spectrometry protocol [48], an overall air flow delivery variability less than 1$\%$, and a rise time (10$\%$–90$\%$ of the final value) close to 300 ms.

2.4. Visual stimuli

For the visual stimuli of neutral faces, we used the Chicago Face Database [49]. We chose 128 different actors showing a neutral expressions. Particularly, we selected a balanced number of actors across gender, age, and ethnicity, to mitigate potential confounding effects related to the intrinsic characteristics of the actors. The pictures were presented in a completely randomized order with respect to the olfactory condition. The visual stimuli were shown on a 15ʹʹ laptop screen with a refresh rate of 60 Hz and a resolution of $1920\times1080$. The pictures' size was 15 cm in width and 10 cm in height, and were displayed at about 50 cm from the eyes of the participants, resulting in a visual angle of about 17∘.

2.5. Experimental protocol

The experimental protocol was divided into two parts. In the first part of the experiment, we collected the subjective ratings of the odorants. To this aim, we presented the olfactory stimuli (i.e. n-butanol, banana and isovaleric acid) to the volunteers in a randomized order for a duration of 10 s, and we asked them to evaluate them in terms of valence (from −2 to +2) and arousal (from 1 to 5) according to the self-sssessment Manikin (SAM) test [50]. Two participants were not able to perceive any of the administered odorants and have been excluded from the experiment, resulting in a final panel of 20 volunteers (age 26 ± 3, 5 females). In the second part of the experiment, we designed an experimental protocol comprised of 128 trials. As schematically reported in figure 1, each trial consisted of: (1) 3 s of dark gray background; (2) 3 s of a dark gray background with a white fixation cross; (3) 1.5 s of neutral face image presentation; (4) 6 s of inter-trial rest. Within each inter-trial rest, subjects were asked to evaluate the facial expression in terms of valence and arousal according to the SAM test, through an interactive interface. The facial images were presented in combination with clean air or one of the three different odorants rated in the first part of the experiment, for a total of 128/4 = 32 trials for each olfactory condition. Each odorant was delivered starting from the onset of the dark gray background to the end of the visual stimulus. A wash-out with clean air was performed during the inter-trial rest. We presented both visual and olfactory stimuli in a randomized order through the PsychoPy software [50]. A custom-made analogic front-end circuit controlled the precise timing for the administration of both the visual and olfactory stimuli.

Figure 1. Schematic illustration of the experimental protocol. Each trial consisted of: (1) 3 s of dark gray background; (2) 3 s of fixation cross; (3) 1.5 s of neutral face image presentation; (4) 6 s of inter-trial rest where subjects rated the valence and arousal of neutral faces according to the self-assessment manikin (SAM) test. A random odor among banana, n-butanol, isovaleric acid and clean air was delivered starting from (1) to the end of (3). A wash-out with clean air was performed during (4).

Standard image High-resolution image 2.6. EEG and EDA acquisition

EEG signal was acquired using a high-density 128-channel geodesic EEG System 300 from Electrical Geodesic, Inc. (EGI). Electrodes were grounded through two additional channels placed between Cz and Pz and referenced through Cz. We always kept electrode impedances below 20 $\textrm\Omega$ during the acquisitions. EEG was acquired at the sampling frequency of 500 Hz.

EDA was acquired using a Shimmer3 GSR+ unit (Shimmer, USA) at the sampling frequency of 250 Hz. We recorded EDA through a pair of Ag/AgCl electrodes placed on the proximal phalanx of the first and second fingers of the non-dominant hand, respectively.

2.7. EDA-driven sympathetic activity estimation

We implemented a procedure based on the analysis of EDA signal to estimate the occurrence of enhanced sympathetic responses associated with the visual presentation of faces. A well-known problem in EDA analysis concerns the temporal overlapping of consecutive SCRs, which may hamper the association between a given response and its potential triggering stimulus [39, 44, 51]. To address this issue, we adopted the cvxEDA [44] model to recover an estimate of the SMNA from the observed EDA responses. Specifically, cvxEDA considers that each SCR is preceded in time by sparse and discrete bursts of SMNA. These bursts are characterized by a higher temporal resolution compared to phasic activity, and can thus be exploited to identify the time instants at which peripheral sympathetic responses evoked by faces occur [44]. Accordingly, we assumed that visual stimuli eliciting an enhanced peripheral sympathetic response could be followed by the occurrence of an SCR and, thus, a non-zero SMNA neural burst.

Operationally, for each subject and for each odor condition (i.e. pleasant, unpleasant, neutral, air), we undersampled the EDA to the sampling frequency of 50 Hz, and we performed a Z-scoring on the data [44]. Then, we applied cvxEDA to obtain an estimate of the SMNA. We set the sparsity parameter of the model to $8 \times 10^$, as a trade-off between noise-suppression and distortion of the solution. Indeed, larger values of this parameter yield to a higher sparsity of SMNA responses and thus to a stronger suppression of noise-induced spikes, but also more attenuation of true physiological responses. On the other hand, smaller values yield to a less distorted but noisier solution [44]. We set the other model's parameters at their default value. We identified discrete events of enhanced sympathetic arousal associated with the presentation of visual stimuli as those epochs having non-zero SMNA bursts occurring in the (1–5) s after stimulus onset. This choice is supported by several studies indicating that a stimulus-evoked SCR is observed to occur within that range of latency after stimulus' onset [39, 52]. A schematic illustration of the procedure is reported in figure 2.

Figure 2. Schematic illustration of the proposed procedure to assess stimuli-related peripheral sympathetic responses from the electrodermal activity (EDA) signal. The raw EDA was preprocessed by downsampling it at the frequency of 50 Hz and applying a Z-score transformation in order to be given as input of the convex-optimization-based EDA (cvxEDA) algorithm. Sudomotor nerve activity (SMNA) estimates are then epoched in the (0–8) s interval with respect to the onset of each visual face stimulus (i.e. 0 s). Significant sympathetic responses to the stimuli are identified as non-zero SMNA bursts occurring in the (1–5) s interval after stimulus' onset.

Standard image High-resolution image

For each odor condition, we then extracted: (1) the number of epochs with/without a stimulus-related sympathetic response (nSymp), (2) the average latency of sympathetic responses, (3) the average amplitude of the SMNA, computed as the average amplitude across epochs and then over time, and (4) the subject-average SCR generated by the SMNA bursts.

2.8. EEG preprocessing

We preprocessed the EEG signal using EEGLAB [53]. First, we filtered the data with a zero-phase low-pass antialiasing filter and then undersampled it to the sampling frequency of 100 Hz. Afterwards, we applied a zero-phase high-pass filter at the cutoff frequency of 0.1 Hz to improve data stationarity. We removed flat and poorly correlated channels by exploiting the method presented in [54]. Specifically, each channel was compared with its reconstructed version obtained from the spherical interpolation of its neighbors and was removed if the correlation coefficient was less than a user-defined threshold. Here, we used a correlation threshold of 0.8. After visual inspection, we recovered the removed channels through spherical interpolation, and we re-referenced the data to its average. We decomposed EEG data through independent component analysis (ICA) [55], and we removed independent components (ICs) resembling artifact activity (e.g. muscular, eye-blinks and other sources of noise) through visual inspection of their associated time course, scalp map, and power spectrum. For each subject, we then epoched the EEG signal in the (−200, 1000) ms range with respect to the onset of visual stimuli (i.e. 0 ms), and we visually inspected the data to remove epochs containing residual artifact activity. An average number of 125 epochs (minimum retained: 99; maximum retained: 128) was retained over the subjects. Finally, clean EEG epochs were corrected for their baseline (i.e. from −200 ms to 0 ms) and low-pass filtered at the cutoff frequency of 30 Hz using default settings in EEGLAB. We grouped the epochs according to the odor condition, and we further distinguished them based on the presence or absence of a peripheral sympathetic activation (hereinafter referred to as Symp condition). Accordingly, we obtained subject-average ERPs for a total of 8 conditions: clean air, n-butanol, banana, and isovaleric acid, with/without the presence of a sympathetic response.

We focused on relevant ERP components associated with the visual processing of faces, as well as their typical regions of interest (ROIs) on the scalp, based on the previous literature [8, 10, 11, 15, 16, 56, 57] (see table 1 for details). Specifically, we extracted: (1) the P1 component in the 120–160 ms interval after stimulus onset in the occipital region around O1 and O2, (2) the N1 component at the same latency as P1, in the central region around Cz, (3) the right/left N170 component in the 160–200 ms interval after stimulus onset in the parieto-temporal regions near P7 and P8, respectively, and (4) the Vertex Positive Potential (VPP) in the 160–200 ms interval after stimulus onset in the same ROI of the N1 component. Moreover, for a comparison with the previous literature, we also identified: (5) the P2 component at 220–300 ms after stimulus onset in the same ROI of P1, and (6) the late positive potential (LPP), a sustained positivity at 300–500 ms after stimulus onset around Pz. For each subject and for each of these components, we extracted the mean amplitudes across the respective time intervals and ROIs.

Table 1. Channels region of interest (ROI) and time interval (ms) for each of the investigated event-related potential (ERP) components (i.e. P1, N1, right N170 (rN170), left N170 (lN170), vertex positive potential (VPP), P2, and late positive potential (LPP)). The amplitude of each ERP component was obtained as the average over the corresponding ROI and time interval reported in the table. Channels' name is reported according to the geodesic EGI 128-channels cap.

 ChannelsTime intervalP1E59 E65 E66 E70 E71 E76 E83 E84 E90 E91120–160 msN1E6 E7 E13 E30 E31 E37 E54 E55 E79 E80 E87 E105 E106 E112 E129120–160 msrN170E96 E97 E101 E102 E108160–200 mslN170E45 E46 E50 E51 E58160–200 msVPPE7 E13 E30 E31 E37 E54 E55 E79 E80 E87 E105 E106 E112 E129160–200 msP2E59 E65 E66 E70 E71 E76 E83 E84 E90 E91220–300 msLPPE52 E53 E60 E61 E62 E67 E77 E78 E85 E86 E92300–500 ms2.9. Statistical analysis on self-report, EDA and ERP data

We investigated for possible differences in the valence and arousal subjective ratings of olfactory stimuli through multiple Wilcoxon sign-rank tests, with a significance level of α = 0.05. P-values were adjusted for multiple comparisons testing with Bonferroni correction. We further tested for the influence of the odors and sympathetic responses on the valence and arousal ratings of faces using linear mixed effect models (LMMs). Given the lack of consensus regarding the most appropriate use of LMMs in the literature (see e.g. [58, 59]), we built our models according to the following rationale: first, we included all of the fixed effects that allowed the model to converge, leading to a full-factorial model. Concerning random effects, we only included those that presented a correlation $|\rho|\lt 0.80$ with other random effects. This, one the one hand, avoided problems associated with model converge and overfitting. On the other hand, it allowed to obtain reliable estimates from the analysis of fixed effects [60]. The significance of each effect was estimated using the Satterthwaite approximation for degrees of freedom. Accordingly, we built two distinct random-intercept LMMs, modeling the valence and arousal ratings of faces, respectively, with Odor (i.e. clean air, butanol, banana, isovaleric acid), Symp (i.e. presence/absence of a sympathetic response) and their two-way interaction as fixed effects (α < 0.05.):

Equation (1)

where $\textrm_i$ is either the valence or arousal subjective rating of faces, $\textrm*\textrm$ is the two-way interaction between fixed effects, and $(1|\textrm)$ is the random intercept for each subject. We could not include any of Odor and Symp among random effects, as they resulted in over-threshold (i.e. $|\rho|\gt 0.80$) multiple collinearities that hampered model convergence.

Concerning the physiological data, we tested for a significant effect of the odors on the SMNA. To this aim, we conducted a $1\times4$ within-subject repeated-measure ANOVA on both the amplitude and the latency of SMNA responses, with Odor as the main effect (α = 0.05). Multiple comparison testing was controlled with false-discovery rate (FDR) for multiple testing under dependency [61]. We further tested for an interaction between the occurrence of sympathetic responses and the olfactory stimuli through a within-subject two-way ANOVA on the number of EDA epochs grouped by odor condition (i.e. clean-air, n-butanol, banana, isovaleric acid) and the presence/absence of an SMNA response, respectively (α = 0.05).

Finally, we tested for a significant effect of contextual odors and sympathetic responses on the amplitude of the ERP components described in section 2.8. To this aim, we followed the methodology explained previously to model the amplitude of each ERP component through a random-intercept LMM, with Odor (i.e. clean air, butanol, banana, isovaleric acid), Symp (i.e. presence/absence of a sympathetic response) and their two-way interaction as fixed effects:

Equation (2)

Where $\textrm_i$ is the amplitude of the ith ERP component, $\textrm*\textrm$ represents the two-way interaction between fixed effects, and $(1|\textrm)$ represents the random intercept for each subject. We could not include either the effect of Odor or the effect of Symp among random effects, as they showed over-threshold (i.e. $|\rho| \gt 0.80$) multiple collinearities. Indeed, on the one hand, that would not have allowed the model to converge. On the other hand, collinearities would have compromised the results of fixed-effects analysis. Post-hoc comparisons were conducted with pairwise t-tests, and the resulting p-values were adjusted with the Bonferroni correction (α = 0.05).

2.10. Effective connectivity analysis with DCM

The analysis of effective connectivity was carried out using SPM12 [62].

We used DCM for ERPs [22, 63] to investigate the modulatory effect of odors and sympathetic responses on the effective connectivity. Specifically, DCM models the observed ERPs by combining a physiologically plausible neuronal model of interacting cortical regions, and a spatial forward model that maps the cortical activity to observed EEG data. Each region is described by three interconnected neuronal subpopulations: interneurons, spiny stellate cells, and pyramidal neurons. Regions are coupled to each other through extrinsic connections, that are distinguished into forward, backward, and lateral according to the hierarchical organization of the cortex [64]. The effect of administered sensory stimuli on neuronal dynamics is accounted for through specific input connections modeling the afferent activity relayed by subcortical structures to the spiny stellate layer of target cortical regions [22, 63]. Notably, such inputs are the same for each experimental condition. Accordingly, differences among ERPs due to either contextual or stimulus attributes are explained by modulatory coupling gains on the connection strengths [22, 63]. The activity of pyramidal neurons from each region is then projected to the EEG channels through an electromagnetic forward model which accounts for the field spread effects in the head.

2.10.1. Network specification

The DCM for ERP framework explains ERP dynamics as arising from a small number of brain regions [65]. The selection of which brain regions to include in the network for DCM analysis can be made using either prior knowledge from the literature or source reconstruction techniques. Here, we adopted a group source reconstruction approach based on multiple sparse priors (MSP) implemented in SPM12 [66]. Particularly, while MSP has been shown to potentially reduce localization error with high-density EEG systems [66], group-inversion yields better results compared to individual inversions by introducing additional constraints on the set of sources explaining subjects' data [67]. Operationally, we used channels' position co-registered to the MNI standard head template as provided by SPM. Then, we inverted ERPs activity on a cortical mesh comprising 8196 dipoles in the time range from −200 to 400 ms. For each subject, we then created contrasts of log power differences in the 0–400 ms time range, collapsed over the experimental conditions, against the prestimulus window (i.e. −200 to 0 ms). These contrast images were smoothed with an 8 mm Gaussian kernel to create a 3D volumetric NIFTI image. Images were entered into a one-sample t-test design in SPM12 and we tested for significant changes in power with respect to the prestimulus through an F-contrast. Significant regions ($p \lt $ 0.05; family-wise error rate corrected at the cluster level) were labeled according with the automated anatomical labeling (AAL) atlas [68].

2.10.2. DCM subject-level connectivity analysis

We performed DCM analysis on the subject-average ERPs relative to the odor conditions, and with and without a sympathetic response. Accordingly, for each subject, we specified a factor analysis on the single-subject connectivity, with Odor and Symp as between-trial factors. We reduced the data to the first four principal components (PC) or modes of EEG channels' mixture. Such a choice is a trade-off between the computational cost of DCM model inversion and the percentage of variance explained by the data [6971]. Notably, reducing ERP data to their first four PC is indicated as sufficient to capture the components of interest [63]. We focused model inversion on the 0–400 ms time window with respect to the stimulus onset, through a Hanning window. Finally, we adopted the ERP neural mass model (NMM) to model temporal dynamics within and between the network sources.

Concerning the forward model, we modeled the spatial activity of brain sources as equivalent current dipoles (ECD option in SPM12) on the cortex. The passive volume conduction effects on the dipoles' electric field were modeled through a BEM model of the head made of three layers, i.e. cortex, skull and brain, whose conductances were set to 0.33, 0.0042 and 0.33 S m−1, respectively.

We allowed the effect of both Odor and Symp to modulate all the connections of the network, including self-inhibitory effects on each node.

2.10.3. PEB group-level connectivity analysis

We inferred th

留言 (0)

沒有登入
gif