Development of the paternal brain in expectant fathers during early pregnancy

The human parenting brain network mediates caregiving behaviors. When exposed to the stimuli of their infants, compared with non-parents, both fathers and mothers exhibit distinct patterns of neural activation. As human males, relative to females, do not undergo robust physiological changes during pregnancy, when and how the paternal brain networks begin to form remains unclear. Thus, using functional MRI, we examined brain activation in response to infant-interaction videos in two groups, childless males and first-time expectant fathers during their partners' early pregnancy before remarkable changes in their partners' appearances commenced. Multivoxel pattern analysis revealed that expectant fathers' left anterior insula and inferior frontal gyrus showed incipient changes in response to parenthood during early pregnancy. Furthermore, these changes were associated with several paternal traits, such as a negative image toward parenting. Such external factors might influence the paternal brain's development during early pregnancy.


Introduction
Decades of research have established the primacy of maternal care as the vital developmental determinant that predicts a broad range of outcomes across adult lifespan, including physical and psychological health problems, as well as cognitive and social development ( Belsky and Fearon, 2002 ;Feldman, 2015 ;NICHD, 2006 ). Due to the increased numbers of fathers involved in direct caregiving due to contemporary sociocultural changes, scientific interest in both human and animal paternal caregiving has been increasing. Researchers have investigated the neural mechanisms and effects of paternal care, measuring brain neurochemistry, genetic and optogenetic targeting of neural circuits to alter animal paternal behavior, neural networks shaped by parental caregiving, and so on ( Dulac et al., 2014 ;Bales and Saltzman, 2016 ;Rilling and Young, 2014 ). Human brain imaging studies have revealed that human males without children show no preferential brain activation when exposed to infant-related (e.g., infant cry) and non-infant-related stimuli (e.g., pure tones), whereas fathers do exhibit activation of several brain regions in response to the stimuli of their infants ( Seifritz et al., 2003 ;Kim et al., 2014 ). Such parent-specific brain responses have been referred to as the global human caregiving network ( Abraham et al., 2014 ;Feldman et al., 2019 ) or the parental brain ( Abraham et al., 2016 ). This network in- ( Gordon et al., 2010 ), whereas a lower testosterone level is associated with more optimal father-infant social interaction ( Weisman et al., 2014 ). These findings raise further important questions regarding how and when such paternal phenotype variability emerges in human males.
Regarding these questions, several longitudinal studies have suggested that some males may already start to acquire paternal phenotypes before childbirth. Expectant fathers showed decreased testosterone levels during the first to third trimesters of their partners' pregnancy ( Edelstein et al., 2017 ;Saxbe et al., 2017 ); however, three major patterns of testosterone change from prepartum to postpartum emerged, with great individual variation ( Berg and Wynne-Edwards, 2001 ). Importantly, individual differences in fathers' oxytocin and testosterone levels after childbirth are associated with the brain's response to infant stimuli ( Abraham et al., 2014 ;Kuo et al., 2012 ). These differences imply that hormonal dynamics, such as modulation of oxytocin and testosterone levels before childbirth, may influence paternal brain behavior development and variability in fathers starting to experience caregiving.
Given these possibilities, we should carefully reconsider previous studies revealing the formation of parental phenotypes in human males. Focusing on the effect of postpartum caregiving experiences, most studies have tried to identify distinct functional, physiological, and anatomical differences by comparing fathers and childless males ( Seifritz et al., 2003 ;Kim et al., 2014 ;Abraham et al., 2014 ). However, if there is considerable variability in the dynamic acquisition of paternal phenotypes among individuals, irrespective of actual experiences of caregiving, or even before childbirth, a simple group comparison between fathers and males not planning to have infants would lead to overlooking the processes. The time at which first-time fathers' parental brain networks begin to develop, and when variability can be observed, is still unclear.
Therefore, in this study, we investigated the shaping processes of the human male parental brain network, especially during their partners' early stage of pregnancy. These fathers had no direct caregiving experience (i.e., first-time expectant fathers), and their partners' bodily appearances had not yet changed remarkably. First, at the group level, we compared brain activities toward the visual stimuli of infantcaregiving behaviors between expectant fathers and males without plans to have infants (control) to investigate to what extent expectant fathers' parental networks would differ from those of the control group. Second, we examined the hormonal and behavioral profiles of expectant fathers during early pregnancy to identify possible factors related to variability in the acquisition of the parental brain network. We measured two hormones, testosterone and oxytocin. As previously mentioned, testosterone profiles of human men with pregnant partners show time-dependent changes ( Saxbe et al., 2017 ) and variability among individuals ( Berg and Wynne-Edwards, 2001 ). Furthermore, the paternal oxytocin level is thought to be associated with the parental brain network and parenting behaviors ( Abraham et al., 2014 ).
We showed that during early pregnancy, expectant fathers' parenting brain networks are still generally similar to those in childless males; however, some specific neural responses could be identified with impending parenthood.

Ethical considerations
This study's experimental protocol and procedures were approved by the Ethics Committee of Kyoto University Unit for Advanced Study of Mind (30-P-8). All participants provided written informed consent before participating in the study.

Participants
Including participants from all socioeconomic statuses, 36 human males with no children (control group, CG; mean age = 31.1 ± 6.2 years) and 36 first-time expectant fathers ( "papa " group, PG; mean age = 33.1 ± 5.7 years; mean partner's gestational age (GA) = 20.9 ± 6.2 weeks) were recruited. All men from the PG and six from the CG were married or living with their romantic partners. All subjects were of East Asian ethnicity (Japanese, n = 70; Korean, n = 2), had normal or corrected-to-normal vision, and were all right-handed. No participant reported any major medical illnesses, major psychiatric disorders, or neurological illnesses. Data for one subject in the PG were excluded from analyses because the partner was in the late stage of pregnancy (GA = 39 weeks), and data for one subject in the CG were excluded because of an inability to complete the fMRI recording session ( n = 35 for each group, see Table 1 ). Data for two subjects in the PG were further excluded from the analyses because of excessive movement artifacts during fMRI recordings. Therefore, the final sample consisted of 33 men in the PG and 35 in the CG.

Behavioral characteristics
Through a questionnaire, we collected the following demographic data: (a) socioeconomic status, including household income, educational background and history, number of family members in the household, and (b) weekly work/study time.
In addition, through standardized questionnaires, we obtained the following participant psychological characteristics: (c) State-Trait Anxiety Index using the State-Trait Anxiety Inventory (STAI-trait) ( Spielberg, 2010 ), and (d) depression index using the Beck's Depression Inventory (BDI) ( Zich et al., 1990 ). These two psychological traits have been implicated in the heterogeneity of parenting behaviors in both mothers and fathers ( Wee et al., 2011 ), as shown by a decrease in positive parenting style and an increase in negative parenting style ( Wilson and Durbin, 2010 ).
We also collected the following behavioral data related to family and parenting: (e) positive and (f) negative image toward parenting ( Inoki and Minori, 2011 ), including 12 items for a positive image (e.g., parenting involves pleasure) and 15 items for a negative image (e.g., parenting is hard), with each item scored on a 5-point Likert scale; (g) partner relationships (Dyadic Adjustment Scale, DAS, ( Spanier, 1988 ) which was only answered by married participants and those living with their romantic partner) with 32 items for assessing marital relationship quality; (h) fetal-paternal attachments (The Paternal Antenatal Attachment Scale, PAAS, ( Condon, 1993 ) only for PG subjects) measuring strength of attachment with 19 items, each scored on a 5-point Likert scale. Scales from (e) to (h) were scored consistently with previous studies. Ceiling/floor effects were observed for several items on the (e), (f), and (h) scales, defined by whether the mean score of ± 1 standard deviation exceeded the range of possible values. Items indicating a ceiling/floor effect were excluded, and mean scores were used for analysis.
Since parenting experiences influence paternal brain activity after birth ( Abraham et al., 2014 ;Feldman et al., 2019 ), earlier parenting experiences might affect activation of the male parental brain. Thus, interviews were conducted to determine whether participants (both PG and CG) had any earlier parenting experiences, such as: (i) past experience of caregiving -whether participants had taken care of babies or young children before adulthood (until the age of 20) -scored 1 or 0; (j) past experience of caring for pets -whether participants had taken care of pets before adulthood (until the age of 20), scored 1 or 0; and (k) recent experience of interaction -whether participants had interacted with infants or young children within the last 2 years -(1 or 0). During interviews and with permission, the participants' voices were recorded. After each interview, one experimenter coded the recorded data offline. Of the data, 25% were also coded by another adult who did not know the study's purpose. Intercoder reliability was found acceptable for these three items: (i) = 0.89, (j) = 0.81, and (k) = 1.00.
The behavioral data is summarized in Table 1 .

Hormone data collection and analysis
Prior to saliva collection, subjects rinsed their mouths, and 10 min later, chewed on an oral cotton swab (Salimetrics, State College, PA, USA) placed sublingually for three min. This process was repeated to obtain a backup sample. The first samples were used for the calculation of each individual's hormonal levels. However, if the first sample contained insufficient saliva or if the first sample's hormone levels were outliers, the second sample's value was used as the individual value. Samples were stored at − 80 °C until assayed. Salivary oxytocin was assayed with a commercial kit (ADI-900-153A-0001; ENZO Co. Ltd., Tokyo, Japan), following the manufacturer's protocol. In brief, 240 μL of saliva was dried using a Speedvac evaporator at room temperature for 3 h and reconstituted in 240 μL of assay buffer, of which 100 μL was used for the assay. Salivary testosterone was measured by ELISA following our previous reports ( Shimozuru et al., 2008 ;Aoki et al., 2010 ); 25 μL of saliva was used in the assay using testosterone-3-CMO -HRP (FKA101; COSMO Bio, Tokyo, Japan) and specific anti-testosterone serum (FKA102-E; COSMO Bio). All intra-and inter-assay coefficients of variation were lower than 15%.
We could not detect oxytocin from the samples of two subjects from the PG and two subjects from the CG. Thus, for the analysis of oxytocin, these subjects were not taken into account.

Task protocol
While inside the MRI scanner, subjects were presented with silent videos of a male model performing actions from a first-person perspective ( Fig 1 ). Subjects were asked to attend to the model's hand movements. The videos were two infant-interaction videos (S1: playing with an infant; S2: changing diapers), and two control videos not showing infant interaction (C1: opening a box and removing a tripod from it; C2: wrapping a box with plastic). Control videos were selected roughly to match movements in the corresponding infant-interaction videos. In the infant-interaction videos, the frame was set to exclude the infant's eyes as previous studies have shown that infant facial features might evoke different responses in human males (i.e., own-vs. other-child ( Kuo et al., 2012 ;Mascaro et al., 2013 ); because of increased response due to facial resemblance ( Platek et al., 2004 ;Platek et al., 2005 )). Therefore, this decision was made to minimize confounding effects in our subjects. The control videos were then framed accordingly.
After an initial 24 s of rest and signal stabilization, each video was presented once per run for 30 s, followed by 30 s of rest. The order of presentation was pseudo-randomized, with the caveat that the infant interaction and its respective control video were always presented consecutively (for example, one run would be S1-C1-S2-C2, while the other was S2-C2-C1-S1). Each subject completed two runs of this task. As part of a different experiment, subjects performed an auditory task immediately after the two video runs. Data from this auditory task were not included in this report.
High-resolution structural T1-weighted images were collected using a three-dimensional magnetization-prepared rapid acquisition with gradient echo sequence, with the following parameters: repetition time = 2250 ms, echo time = 3.51 ms, field of view = 256 mm, matrix size = 256 × 256, flip angle = 90°, voxel size = 1 × 1 × 1 mm, slice gap = 0; number of slices = 208 axial slices; slice order = interleaved. fMRI data were preprocessed using MATLAB R2018b (Math-Works, Natick, USA) and the SPM software package (SPM12 v7487, https://www.fil.ion.ucl.ac.uk/spm/ ). Each run's first eight volumes were discarded to allow for signal stabilization. Functional images were corrected for acquisition time, resliced, realigned, and normalized to match the MNI template brain and were spatially smoothed with a Gaussian kernel of full width at half maximum (FWHM) of 4 mm. Subjects with movement artifacts greater than 3 mm within the runs were discarded.
We modeled each subject's response to the stimulus using a general linear model (GLM), in which stimulus blocks were defined as predictors and convolved with the standardized model of hemodynamic response function. For analysis, we focused on three contrasts: S1-C1, S2-C2, Fig. 1. Example of a typical run inside the MRI scanner. Participants were asked to view alternating silent infant interaction (Stim) and matched control videos. Videos were shown for 30 s, with 30 s of rest afterwards. The presentation order was pseudo-random. and the overall S-C (combining S1 and S2 and combining C1 and C2). Contrasts were defined as stimuli greater than the control condition.

fMRI data analysis 2.7.1. Whole-brain analysis
Whole-brain analysis was conducted using SPM and a second-level design. Individual contrast data for S-C for all subjects were included in a one-sample upper-tailed t -test model to test the null hypothesis that the group-level observed signal in a given voxel cluster would be due to a non-task effect. Significance was assessed by initially setting a voxel-wise uncorrected significance level of p < 0.001, with a minimal cluster size of 20 voxels, and then selecting suprathresholded clusters with a family-wise error-corrected significance threshold of p < 0.05. Using the WFU Pick Atlas toolbox for MATLAB ( Maldjian et al., 2003 ;Maldjian et al., 2004 ) and the AAL atlas of the human brain ( Tzourio-Mazoyer et al., 2002 ), we matched surviving clusters to known anatomical areas (Supplementary Table S1).

Region of interest (ROI) analysis
To study the functional response to infant-interaction stimuli in greater detail, we conducted ROI analysis. To reduce selection bias, we created functional ROIs based on previous research using a similar stimulation paradigm ( Abraham et al., 2014 ), corresponding to parenting network areas. Most suprathresholded clusters of activation in the whole-brain analysis match the ROIs described by Abraham et al. These ROIs are the ventro-anterior cingulate cortex (vACC), amygdala (AMY), left inferior frontal gyrus/anterior insula (IFG_L), left middle insula (INS_L), bilateral lateral frontopolar cortex (LFPC), STS, temporal pole (TP), and ventromedial (vm) prefrontal cortex (PFC). Because Abraham et al. did not find significant activation in the right inferior frontal gyrus and right insula, they did not include those two areas' coordinates; however, we did find a significant cluster matching those areas. Thus, we took coordinates for these two areas (IFG_R and INS_R) from Mascaro et al. ( Mascaro et al., 2014 ) because, in their publication, Mascaro and colleagues describe men's functional response to infant cry stimuli. Although their paradigm was different, the lack of research on this topic left us with no suitable alternative.
We created 15 mm-side box-shaped functional ROIs (fROIs) around the peak coordinates for each of the areas mentioned above. Previous research showed no differences between the left and right sides of the AMY, STS, vACC, vmPFC, LFPC, and TP. Thus, we joined the lateral sides of these fROIs into a single one, and fROIs were masked to feature only gray matter.
From each of these fROIs, we extracted blood-oxygen-leveldependent (BOLD) signals at the individual voxel level for each subject and for each contrast (S1-C1, S2-C2, and S-C) using the Brain-DecoderToolbox2 software package ( https://github.com/KamitaniLab/ BrainDecoderToolbox2 ) for MATLAB. Preprocessed data were shifted by two volumes (6 s) to account for hemodynamic delay, adjusted for movement parameters, and linearly de-trended. We then calculated the mean percentage of change from baseline across all voxels in a given fROI as the metric for activation of that fROI. Baseline was defined as the average BOLD signal during rest periods across the experiment. Finally, we removed outliers at the voxel level (threshold of three times the standard deviation).
To test whether each fROI contrast correlated with any behavioral traits or hormonal levels, we used a linear regression model for each fROI-covariate pair. We then evaluated regression coefficients using an F-test with the null hypothesis that the model's fixed-effects coefficients equal 0, correcting p-values by the number of fROIs.

MVPA
MVPA was implemented using libsvm ( Chang and Lin, 2011 ). We trained a support-vector machine (SVM) classifier with a linear kernel, using each subject's non-averaged (across the fROI) and non-spatially smoothed voxel data from the S-C contrast from each fROI to perform binary classification of a target label (see below). To validate the model, we performed leave-one-out cross-validation. In short, for each fROI, we trained a model using voxel data of a given fROI from all subjects but one and used the excluded subject as test data. This procedure was repeated for each subject until all subjects had served once as test data. Then, we averaged the prediction accuracy for each of the test data to obtain the classifier's mean prediction accuracy for that fROI. We also show the 95% confidence interval, calculated using the normal approximation method, for the binomial confidence interval. High accuracy implied that the model could separate multivoxel data between subjects who showed a certain paternal characteristic and those who did not and, thus, could establish a relationship between neural activation patterns and paternal traits. Significance was assessed with an upper-tailed binomial test for the null hypothesis that the resulting classification accuracy emerged from a binomial distribution with a mean equal to the chance level (50%), corrected for the number of fROIs. Additionally, we calculated the true positive rate (TPR) and true negative rate (TNR) for each ROI -target label pair, defined as the ratio between the number of test labels that were correctly classified over the total number of labels that were classified as such (in our case, TPR refers to classification of low behavioral characteristics, whereas TNR refers to high behavioral characteristics; see next paragraph).
The predicted target labels were of two kinds: paternity status (same as group, either PG or CG) and behavioral characteristics obtained through surveys and interviews (see above) and hormonal levels. Prediction of the former was straightforward because paternity status was already a binary variable. However, most of the other traits were not. Because SVM works best with binary classification, we binarized all behavioral traits and hormonal measurements by calculating the median of each across all subjects and assigning the descriptors low (0) to those values under the median and high (1) to those higher or equal to the median. We used bootstrapping to deal with imbalances in sample sizes, if any.

Voxel-based morphometry (VBM) analysis
VBM analysis was conducted using CAT12 v12.6 ( http://dbm.neuro. uni-jena.de/cat12/ ) and SPM12. This analysis allowed us to study structural changes in the gray matter ( Wright et al., 1995 ;Ashburner and Friston, 2000 ). Each subject's T1 anatomy data were normalized and segmented into gray matter, white matter, and cerebrospinal fluid, following the standard CAT12 pipeline. The total intracranial volume (TIV) for each subject was also calculated. Gray matter was then spatially smoothed with a Gaussian kernel with an FWHM of 4 mm. Whole-brain analyses were conducted using second-level designs in SPM, similar to fMRI data. Individual gray matter data for each PG subject were input into a GLM with GA as the covariate and TIV as the nuisance variable. Significance was assessed by initially setting a voxel-wise uncorrected significance level of p < 0.005 and then selecting suprathresholded clusters with a false discovery rate-corrected significance threshold of p < 0.05. These parameters were set in this manner to mimic the procedure outlined by Kim et al., 2014( Kim et al., 2014, who similarly studied structural changes in fathers' brains during the first months postpartum.

Whole-brain analysis
Using functional MRI, we recorded the neural activations of firsttime expectant fathers (PG) and control childless males (CG) as they watched videos of a man performing actions recorded from the firstperson perspective ( Fig 1 ). Stimuli consisted of four scenes, two showing male-infant social interactions (S1: playing with an infant; S2: changing diapers), and two showing control non-infant-interactions (C1: opening a box and removing a tripod; C2: wrapping a box with plastic). Each control scene's motor movements roughly matched those of the corresponding infant-interaction scene. To identify neural activation patterns toward the visual stimuli of infants common to all participants, we performed whole-brain analyses for infant videos vs. control videos (S-C) for all data, regardless of paternity status (i.e., both PG and CG).
Contrary to previous findings that childless males exhibited little neural response to infant stimuli ( Seifritz et al., 2003 ), we revealed extensive activation patterns common to all participants, consisting of 11 distinct voxel clusters throughout the brain ( Fig 2 , Table S1). These clusters corresponded to the following three axes linked to human parental brain networks: three clusters in embodied simulation (supramarginal gyrus, inferior frontal gyrus, and anterior cingulate cortex), three in emotional processing (ventrolateral and dorsolateral prefrontal cortices and the insula), and three in mentalizing (ventromedial and dorsomedial prefrontal cortices and the STS). The remaining two clusters were not associated with the human caregiving network: one in the occipital lobe (Cluster 2), corresponding to the visual cortex, and the other in the left postcentral gyrus-superior parietal lobule (Cluster 6).
Next, we examined differences between expectant fathers and control male neural activation patterns. Second-level ANOVA analysis re-vealed no significant differences in whole-brain activations between the two groups for the S-C contrast. There were no significant differences between the two groups for each contrast between S1-C1 and S2-C2. These results suggest that expectant fathers do not yet show specific neural responses to infant stimuli, at least during early pregnancy (20.9 weeks GA on average). Table 1 summarizes behavioral traits and salivary hormonal levels of expectant fathers and control males. The only significant betweengroups differences identified by a two-sample two-tailed t -test were in work time (daily work hours: 9.54 ± 2.03 in the PG [ "papa " group] vs. 6.11 ± 3.18 in the CG, p < 0.0001; weekly workdays: 5.29 ± 0.62 vs. 3.77 ± 1.96, p < 0.0001; weekly work hours: 50.54 ± 11.95 vs. 28.49 ± 16.32, p < 0.0001), and household income (¥6-7 M ± 3 M, vs. ¥3-4 M ± 2 M, p < 0.0001). Pearson's linear correlation test across covariates revealed a significant positive correlation between household income and weekly work time ( R = 0.5647, p < 0.05, corrected for number of comparisons) and between trait anxiety and the Beck Depression Inventory score ( R = 0.7147, p < 0.05, corrected for number of comparisons). A significant negative correlation was found between a positive image toward parenting and fetal-parental attachment in the PG ( R = − 0.6756, p < 0.05, corrected for number of comparisons). Despite these differences, no behavioral traits or hormonal levels showed any significant correlation with neural response for the S-C contrast.

ROI analysis
However, whole-brain analysis could not detect smaller or weaker effects in certain areas of the brain, particularly in those belonging to the human parental caregiving network. Therefore, focusing on areas previously identified as part of this network ( Abraham et al., 2014 ;Feldman et al., 2019 ;Mascaro et al., 2014 ), we applied ROI analysis. Functional ROIs (fROIs) were selected to reflect the two main components of the human caregiving network: the mentalizing network and the emotional processing network. Based on these two networks, we selected the following fROIs: bilateral ventro-anterior cingulate cortex, bilateral amygdala, left and right inferior frontal gyrus, left and right insula, bilateral LFPC, bilateral STS, bilateral temporal poles, and bilateral ventromedial prefrontal cortex. From each fROI, we extracted the percentage BOLD signal change over the baseline for the infant videos minus the control videos' contrast (i.e., the average%BOLD signal change {S1 + S2}-the average of%BOLD signal change {C1 + C2}). However, no differences were found between PG and CG in any of the selected ROIs (Supplementary Table S3).
Here, worthy of attention in the PG is that the partner's GA significantly correlated with activation of the left inferior frontal gyrus (linear regression coefficient R 2 = 0.2631, p < 0.05) for the S-C contrast. For each of the S1-C1 and S2-C2 contrasts, we also found significant fits for the linear regression model between GA and the left inferior frontal gyrus (R 2 = 0.3129, p < 0.01) and the amygdala (R 2 = 0.2959, p < 0.05) ( Fig 3 ).

Classification of paternity status
Our findings imply that expectant fathers' parental brain networks might begin to change depending on the experience of their partners' pregnancy. To examine this possibility, we used a more sensitive method, MVPA, to determine whether we could detect some patterns of neural activation characteristics of expectant fathers in the paternal network in response to infant-interaction stimuli. We applied a machine learning method, utilizing an SVM model that performs supervised classification of paternity status using multivoxel data in each fROI in all participants. For each participant and fROI, the algorithm created a model using voxel data in an fROI from all other participants;

Fig. 2.
Whole-brain analysis for the Infant videos > Control videos revealed 11 voxel clusters of activation common to all participants. These clusters corresponded to the following three axes linked to human parental brain networks: embodied simulation (SMG, and IFG), emotional processing (INS), and mentalizing (vmPFC, dmPFC, and the STS). p < 0.001 uncorrected, minimal cluster size k = 20 voxels. Coordinates in MNI space. SMG, supramarginal gyrus; IFG, inferior frontal gyrus; INS, insula; vmPFC, ventromedial prefrontal cortex; dmPFC, dorsomedial prefrontal cortex; STS, superior temporal sulcus; PCC, posterior cingulate cortex.  Asterisks denote a significant difference from the chance level ( p < 0.05, corrected for number of fROIs). (b) Following the same procedure as before, we successfully predicted for the PG the following: binarized negative image toward parenting in the bilateral amygdala (AMY), partner's GA from the bilateral amygdala, weekly worktime from the right insula (INS_R), and recent experience with childcare from the bilateral superior temporal sulcus (STS) and the bilateral temporal pole (TP). From the CG voxel data, we could only achieve a marginally significant greater-than-chance classification accuracy for recent experience with childcare from the TP. Asterisks denote a significant difference above chance level, not between groups (binomial test, * * : corrected p < 0.0001, * : corrected p < 0.05, † : uncorrected p < 0.01). We did not perform statistical comparison between PG and CG groups. then, it would predict the participant's paternity status based on that model. The algorithm was able to classify paternity status from the left insula with relatively high accuracy (mean ± 95% confidence interval: 67.65 ± 11.12, p = 0.029342, TPR = 0.70, TNR = 0.66 , Fig 4 a). With the expectant fathers differing from childless males, this suggests a common neural substrate in encoding patterns for male-infant-interaction stimuli in expectant fathers.

Classification of paternal traits
Previous studies have reported a wide variability in testosterone levels ( Fan et al., 2011 ) and activation profiles in the paternal brain ( Abraham et al., 2014 ). We assumed that some expectant fathers' paternal characteristics, such as behavioral traits and two kinds of hormonal levels (i.e., testosterone and oxytocin), would relate to individual variability in the parenting network's activation. However, we found no significant correlation between any behavioral traits (i. e., parental attitude, anxiety, depression, or fetal attachment) and BOLD signal change in any fROI (see Supplementary Tables S5 -S7). We also found no significant correlation between hormonal levels and BOLD signal change to infant stimuli during early pregnancy (Pearson's correlation, p > 0.05).
Our MVPA findings suggest that rather than the parental network's overall activation, multivoxel patterns provide more information of the paternal brain's early development. Applying an SVM model, we performed supervised classification of binarized behavioral traits (i.e., lowhigh) using multivoxel activation patterns of each of the selected fROIs. In PG, the following characteristics were significantly predicted (above chance level) from different ROIs ( Fig 4 b, Supplementary Table S8): Participants' negative image toward parenting was predicted from the amygdala (mean ± 95% confidence interval: 87.88 ± 11.14, p < 0.0001, TPR = 0.73, TNR = 1.00), partner's GA was predicted from the amygdala (75.76 ± 14.62, p = 0.023, TPR = 0.71, TNR = 0.79), weekly worktime was predicted from the right insula (75.76 ± 14.62, p = 0.023, TPR = 0.60, TNR = 0.89), recent experience with childcare was predicted from the STS (72.73 ± 15.20, marginally significant at p = 0.068, TPR = 0.30, TNR = 0.91), and recent experience with childcare was predicted from the temporal poles (75.76 ± 14.62, p = 0.023, TPR = 0.60, TNR = 0.83). Additional analysis confirmed that, with the exception of the amygdala -GA pair, changes in BOLD signal obtained from univariate analysis did not predict these fathers' behavioral profiles ( Supplementary Fig.  S1). Regarding the ROI -covariate pairs that showed significant effects in PG, MVPA in CG revealed only a marginally significant association between recent experience of childcare with infants and multivoxel patterns in the temporal poles (71.43 ± 14.97, p = 0.083, TPR = 0.88, TNR = 0.58, Supplementary Table S9). These results were consistent with the view that during early pregnancy, expectant fathers show nascent but subtle neural responses to imminent parenthood.

VBM analysis
Previous studies have suggested that pregnancy in human primiparous females is associated with pronounced changes in brain structure, especially extensive gray matter volume reductions across pregnancy ( Hoekzema et al., 2017 ). Based on this, we also conducted a VBM analysis to compare expectant fathers and control male brain structures. However, no differences were found. Next, we applied a GLM to link gray matter volume changes with expectant fathers' gestational age, finding a significant negative correlation between GA and three clusters in the brain, corresponding mainly to the amygdala and hippocampus, thalamus, and right insula ( Fig 5 and Supplementary Table S4). Of note, this finding was consistent with the relationship between fMRI activation patterns in the amygdala and GA, as found by MVPA.

Discussion
In this study, we revealed three main findings. First, human males exhibit preferential neural activation in several areas of the parental caregiving network in response to infant-interaction stimuli when compared with matched control stimuli, regardless of paternity status. Sec-

Fig. 5.
Whole-brain analysis of the relationship between expectant father's gray matter and partner's GA. Using a GLM, we found three clusters with a significant negative effect between gray matter volume and partner's GA. Labeling was performed using AAL atlas of the human brain. AMY: amygdala, HPC: hippocampus, INS: insula, THA: thalamus.
ond, during early pregnancy, expectant fathers' brains begin to show incipient changes in response to imminent parenthood. Third, among expectant fathers, differences in neural activation patterns in areas of the paternal brain may be related to some specific behavioral characteristics, although these differences are not related to their hormonal levels, at least during early pregnancy.
Regarding the first finding, we measured the neural response to infant-interaction stimuli in human males with no children and in expectant fathers during the early stage of their partners' pregnancy. Wholebrain analysis revealed preferential activation in response to infant stimuli in nine different clusters across participants' brains, irrespective of paternity status. All these areas have been previously associated with the human caregiving network ( Abraham et al., 2014 ;Feldman et al., 2019 ). To our knowledge, this is the first study to examine human male brain responses toward infant interaction, including those with no plan to have children. The current study provides evidence that, regardless of paternity status, human males show preferential activation to visual stimulus of father-infant-interaction compared to matched non-infantinteraction videos. It may be possible that some of these activations are due to an object-interaction effect from our tasks; however, given the simplicity of the tasks (extracting an object from a box and wrapping a box with plastic), we do not expect these stimuli to elicit a strong emotional or cognitive response in our subjects, particularly in non-motor areas. Incidentally, we observed some clusters of activation in visual and visuo-motor processing areas in both the infant > control and control > infant stimuli, suggesting that while there are some object-interaction effects, they occur in separate areas.
Next, we found no significant differences in fROI mean activations between expectant fathers during early pregnancy and control males. Expectant fathers' brains were activated in the same manner toward infant-interaction stimuli as those of control males. This supports the view that the paternal brain forms either late in pregnancy or after childbirth when fathers start to accumulate caregiving experiences with their infants ( Feldman et al., 2019 ). In contrast, we found significant associations between the partner's GA and the expectant father's activation of two parenting-related brain areas, the left inferior frontal gyrus, and the amygdala, which are activated in response to infant-interaction scenes, that is, social play with an infant and changing diapers, respectively. These results imply that, at least in some areas, the paternal brain might start to develop during early pregnancy. According to the AAL atlas of the human brain ( Tzourio-Mazoyer et al., 2002 ), the selected left inferior frontal gyrus fROI maps roughly to anatomical areas corresponding to the inferior frontal gyrus with some overlap with the anterior insula; both have been associated with emotional processing ( Fan et al., 2011 ;Craig, 2009). Wittfoth-Schardt et al. ( Wittfoth-Schardt et al., 2012 also showed that fathers exhibit greater activation of the left inferior frontal gyrus/anterior insula when watching images of their children than when watching images of other children. Although not as clearly as after childbirth, the expectant fathers in our study might have started to specifically image imminent caregiving behaviors toward their nascent infants. One explanation for the negative relationship with GA is that as pregnancy progresses, the mental image of fathers of their unborn infants may develop more, in turn leading to stronger differences between this image and the infant in our stimuli, resulting in decreased activation in the inferior frontal gyrus.
The link between GA and amygdala activation in response to the diaper change stimulus may reflect a change in the perceived emotional valence of a core caregiving behavior in expectant fathers. The amygdala has been associated with the response to emotional stimuli of both negative and positive valence; a decrease in the amygdala activation as the pregnancy continues could be attributed to fathers reducing their negativity toward an action that once seemed challenging. It could also indicate a reduction in the positive feelings of expectancy toward changing their own infant's diaper as the reality of their new fatherhood settles in. Another possibility is that, as the pregnancy progresses, an expectant father becomes more accustomed to infants, reducing the emotional salience of our stimuli. Without an individual baseline to compare against, we cannot conclusively determine the mechanism driving this GA-dependent decrease in amygdala reactivity. That said, this finding, paired with the inferior frontal gyrus-GA link described above, provides the first evidence that pregnancy-dependent changes take place in the brains of expectant fathers.
We also successfully separated the two groups based on fine patterns of activity in the left insula, identified using MVPA. This finding suggests that some patterns of neural activation are common to expectant fathers but not to childless males. The left insula has been associated with emotional and cognitive empathy ( Uddin et al., 2017 ) and might mediate some parental behaviors ( Rilling, 2013 ). In addition, compared to childless males, expectant fathers exhibit greater activation in the left and right insula in response to infant cries ( Seifritz et al., 2003 ). Our results with MVPA might signal an incipient change in expectant fathers' brains that will eventually lead to more remarkable differences between the two male groups.
Regarding the third finding, we could correctly classify different paternal traits among expectant fathers using MVPA, which suggests underlying differences in neural activation patterns. Compared to those with partners in mid-pregnancy (GA ≥ 19 weeks), men with partners in earlier pregnancy (GA < 19 weeks) showed different patterns of neural activation in the amygdala in response to our stimuli. Differences were also observed between the amygdala's neural activation patterns common to fathers with a high negative image toward parenting and the patterns common to fathers with a low negative image. Similar differences in activation patterns were also found in the temporal poles and STS (albeit marginally significant) among expectant fathers, whether they had recent experiences with childcare or not, and in the right insula among individuals, depending on their weekly work time. These in-group differences in neural encoding patterns observed only in expectant fathers suggest that the brains of some, but not all, might initiate nascent changes in response to their imminent role as fathers and that such changes might depend on several external psychological and behavioral factors mentioned above. Our findings imply that, from early pregnancy, such external factors might influence the variable acquisition of paternal phenotypes.
The relationship between the recent experience of STS and expectant fathers with childcare should be highlighted. One study reported that parental STS activation in response to own-infant interaction stimuli depends on experience with childcare after childbirth ( Abraham et al., 2014 ). It is possible that the current finding reflects an incipient mechanism that matures after the infant's birth and subsequent parent-child interaction. If so, our findings provide evidence that paternal brain development before childbirth could be modulated by several external factors. Thus, our findings might pave the way for interventions and other aids designed to help men adopt the role of caring and involved fathers.
In contrast, we found no association between expectant fathers' hormone levels (oxytocin and testosterone) and their neural activation during early pregnancy. Testosterone gradually decreases during late pregnancy ( Gordon et al., 2010 ), and oxytocin increases postpartum during interaction with the infant ( Edelstein et al., 2017 ). Therefore, the effect of these hormones on a father's brain development might become apparent from late pregnancy to postpartum, when hormonal levels greatly fluctuate. Further longitudinal studies during the entire pregnancy period are necessary to confirm these assumptions. Additionally, Matsunaga et al. ( Matsunaga et al., 2020 ) showed that baseline oxytocin levels show a high degree of variability. In our current dataset, we also observed high variability in both salivary oxytocin and testosterone levels in our participants. Our sample size may not have been large enough to overcome this variability and reliably detect a relationship between hormone levels and brain activation. A longitudinal study may be able to bypass this issue by measuring changes in hormonal levels across time rather than at a baseline at a single point.
Why do these differences in multivoxel patterns of neural activation occur in expectant fathers but not in childless males? One possibility is that the phenomenon might be related to cortical restructuration, either structural or functional, evoked by the prospect of fatherhood in parenting-related areas. Kim et al. ( Kim et al., 2014 ) reported structural changes in the human male's brain, including decreased gray matter volume in the insula, during the postpartum period (from 2 to 4 to 12-16 weeks after childbirth), suggesting that caregiving experiences can induce cortical plasticity. In contrast, when comparing at pre-and postpregnancy, no significant cortical volume changes occurred in human males ( Hoekzema et al., 2017 ). These two studies suggest that structural changes in fathers' brains begin only after childbirth. However, our structural analysis using VBM revealed small reductions in gray matter in the right insula and amygdala, but not the left insula, as gestation progressed. We assume that a reduction in cortical volume linked with gestation progression might be a precursor to fathers' postpartum cortical restructuration, as shown by Kim et al. ( Kim et al., 2014 ). These structural changes might partially underlie differences between expectant fathers in multivoxel patterns of neural activation in the amygdala in response to infant-interaction stimuli. The differences in multivoxel patterns that we found in other areas might instead reflect functional restructuration. VBM is a technique that measures regional gray matter volume and is sensitive to some, but not all, changes in cortical structure. For example, the rewiring of some neural circuits may not be large enough to produce changes in the local gray matter volume. Thus, MVPA, which reflects functional activation patterns in a group of voxels in response to a stimulus, might be a more effective method for detecting change even in areas where VBM finds no apparent change.
Finally, it is important to discuss the limitations of our study. Our choice of contrasting stimuli, infant-interaction vs. object-interaction, may raise important confounding factors. Primarily because we did not control for social interaction (i.e., infant-interaction vs. adult-interaction), we cannot discard the possibility that the neural activations we observed in the brain areas belonging to the emotional processing and mentalizing networks in response to infant-interaction stimuli are due to a general social effect, rather than a parental-specific effect. A meta-analysis study of neuroimaging studies regarding social perception in humans by Arioli &Canessa, 2019 ( Arioli andCanessa, 2019 ) revealed that areas of the mentalizing network and emotional processing networks are also engaged in a variety of social interaction situations between adults. Thus, the possibility that our whole-brain analysis findings, particularly those of control men, may be due to a regular nonparental social effect cannot be ignored. In contrast, most of the brain areas we reported have also been reported in previous parental brain studies ( Abraham et al., 2014 ;Mascaro et al., 2014 ;Young et al., 2017 ), thus supporting our findings. In addition, our results relating to an association between brain activations, gray matter changes, and GA, or between other parental-related behavioral traits are probably parentalspecific, due to their parental context. Further studies contrasting the differences between the neural responses in men to peer social interaction and infant interaction are needed to shed more light on this matter. Another limitation of our study is our relatively small sample size, particularly with respect to the number of tests conducted. We corrected p values only on the basis of number of ROIs and not according to the number of behavioral and hormonal traits as this would have drastically lowered the statistical power of the present study. Thus, further studies with larger sample sizes are needed to confirm our findings. Moreover, our choice of control group could be considered as an additional limiting factor. Since the control group mostly comprised single men, which was in contrast with the expectant fathers group that had married men, it is possible that the differences observed between the two groups were partially confounded by the discrepancies in the marital statuses of the participants. However, it should be noted that our metric for partner relationship, DAS, showed no significant effect on brain activation in fathers in any of our analyses, suggesting that including married men in the control group would perhaps have not provided any additional benefit. Additional studies exploring the effect of marital status on the development of the paternal brain are needed to address this issue.
In conclusion, our findings support the hypothesis that development of the paternal brain begins during early pregnancy in first-time expectant fathers. Further research is needed to clarify whether these early changes continue linearly throughout pregnancy and after childbirth, and to determine the extent to which individual differences play a role in the paternal brain's development.

Data availability statement
Due to strict privacy concerns, behavioral, hormonal, and MRI data obtained throughout this study can only be made available from F.D.R. upon reasonable request.

Code availability statement
All custom scripts and codes used to generate the results hereby presented in this study are available from https://github.com/fdiazro01/ PapaPro . Main routines used for whole-brain analysis, ROI signal extraction, MVPA, and VBM are available from their respective owners, as described in the Methods section.