1. Introduction. Social cognition, which refers to how individuals process, memorize, and use information in social contexts to explain and predict their own behavior and that of others (Fiske and Taylor (1991)), is a crucial aspect of human functioning and has been extensively studied in the field of psychology and neuroscience (Lieberman (2007), Saxe and Kanwisher (2013)). The use of neuroimaging techniques, particularly functional magnetic resonance imaging (fMRI), has enabled a better understanding of the neural mechanisms underlying social cognition (Saxe and Kanwisher (2013)). Previous studies using fMRI have shown that specific brain regions, such as the medial prefrontal cortex, the temporoparietal junction, and the superior temporal sulcus, are consistently activated during tasks related to social cognition (Castelli et al. (2000), Gallagher and Frith (2003)). While significant progress has been made in uncovering the neural mechanisms underlying social cognition, our understandings of the coordination between brain regions during social cognition and how it relates to individual differences in social behavior remain limited (Adolphs (2009)).

The social cognition study in the Human Connectome Project (HCP) 1 provided a unique opportunity for advancing our understandings of the brain connectivity underlying social cognition. In this study, imaging scans are collected using fMRI from a set of subjects as each subject goes through a sequence of cognitive tasks and rest states. In addition, it also collects subject features such as sex and social covariates (e.g., social distress); see more details in Section 1.1. Based on the imaging scans, a dynamic connectivity network, characterizing activation and deactivation of connections between brain regions during task and rest states, can be constructed for each subject, with nodes corresponding to a common set of brain regions and the edges encoding dynamic functional associations between the regions. From this study it is of fundamental scientific interest to understand which brain regions are co-activated during the cognitive tasks. In addition, it is important to understand whether there are sex differences in brain connectivity during cognitive tasks, and if so, how social covariates influence these differences.

There is some recent literature on modeling a collection of networks, including dynamic networks. However, these methods may not flexibly associate dynamic network connectivity with external covariates while taking into account the structure of the network and smoothness in the dynamic brain connectivity. Specifically, Xu and Hero (2014), Pensky (2019), Zhang and Cao (2017), Zhang, Sun and Li (2020) proposed several approaches based on stochastic block models. These methods cannot associate network connectivity with external covariates. Wang et al. (2017) proposed a Bayesian network model with covariates, which is flexible but can be computationally intensive, especially for large networks or a large number of covariates. Kong et al. (2020), Hu et al. (2021), Zhang, Sun and Li (2023) studied matrix or network response regressions, but they focused on nontime-varying networks. Sun and Li (2017), Hao et al. (2021), Zhou et al. (2021), Tang, Bi and Qu (2020) considered tensor regressions that can be formulated to tackle our problem by stacking the dynamic networks observed at different time points into a tensor, but these approaches could not account for the temporal smoothness in the dynamic brain connectivity.

To model the dynamic brain connectivity in the social cognition study, we propose a new semiparametric dynamic network model for a collection of dynamic networks with subjectlevel covariates. We adopt the form of generalized linear model (GLM) and assume the connectivity between a pair of regions, after a proper transformation, is the sum of two functional components. The first component is the baseline time-varying connectivity shared by all subjects, and the second component involves time-varying slopes and models the effects of subject-level covariates on the time-varying brain connectivity. To estimate the unknown functional coefficients, we consider a nonparameteric estimation via B-spline approximations. Under such approximations we can then write our model in the form of a dynamic network regression, where the response is the dynamic connectivity matrix and the predictors are subject covariates. With the B-spline basis, the baseline connectivity can be characterized using an intercept tensor and the covariate effect using a slope tensor. We assume the intercept tensor is low-rank and the slope tensors are structurally sparse. We discuss the benefit of placing different assumptions on these two tensor coefficients in Section 2.1. These structural hypotheses significantly reduce the number of free parameters, facilitate model interpretability and estimability, and are commonly considered in scientific applications (Bi, Qu and Shen (2018), Zhang, Sun and Li (2023)).

For estimation we propose an efficient alternating gradient descent algorithm with a fast iterative shrinkage-thresholding method to estimate the sparse slope tensor. In Section 3 we demonstrate in simulation studies that our method can accurately estimate the model coefficients and identify nonzero covariate effects, whereas other methods fail to offer accurate estimates. In Section 4 we apply our proposed method to the social cognition study and identify sex differences both in the baseline connectivity and social covariate effects. The majority of our results agree with the existing findings in the neuroscience literature. We also implement an elementwise (i.e., edge-based) method, where the results are highly noisy and lack interpretability, and a method designed for nontime-varying networks (Zhang, Sun and Li (2023)), where the results are highly sparse and cannot identify areas that are known to be engaged in social cognition. Finally, we consider a permutation based procedure to evaluate the identified sex-specific differences from our analysis.

Taken together, our work proposes a new dynamic network regression for analyzing taskevoked brain connectivity with subject-level covariates that exploits the structure in the brain network and the temporal smoothness in the time-varying connectivity. We demonstrate in simulations and real data analysis that the proposed method usually performs better than elementwise methods that model the connectivity between each pair of nodes separately. Next, we discuss in detail the motivating scientific problem and the research questions to be addressed.

1.1. The HCP social cognition study and research questions. The social cognition study in the HCP data collected behavioral and task-related fMRI data from 850 healthy adult subjects. In each session a participate was presented with several short videos of objects (squares, circles, triangles) interacting (Castelli et al. (2000)), and the fMRI data were collected on 274 evenly spaced time points. These videos were developed by either Castelli and colleagues (Castelli et al. (2000)) or Martin and colleagues (Wheatley, Milleville and Martin (2007)). Specifically, two types of video clips were shown to the subjects including mental (objects interact in some way: see Figure S13 in the Supplementary Material, Zhang et al. (2024)) and random (objects move randomly). For each participant, there were five video blocks (three mental and two random), with each video task and rest duration taking up 23 seconds and 15 seconds, respectively. We focus our analysis on the N = 843 subjects who were shown videos in the sequence of mental, mental, random, mental and random. Additionally, social related traits, such as social distress, social support, and companionship, were measured for each subject via self-reported questionnaires; see more details in Section 4.

In our analysis the fMRI data are preprocessed and summarized as a 68 × 274 spatialtemporal matrix for each subject using the Desikan-Killiany Atlas (Desikan et al. (2006)) with n = 68 regions of interest (ROIs, see Table S4). As each subject goes through various tasks and rest states during the scanning session and activation/deactivation of brain regions measured via fMRI are typically lagged (Schölvinck et al. (2010)), it is more appropriate to study the brain connectivity as a dynamic network. Specifically, for each subject the dynamic network is constructed by calculating a sequence of connectivity matrices over T sliding windows, each summarizing the connectivity between 68 brain regions in a given window. While there are many choices of connectivity measures (Smith et al. (2013)), the most commonly used one is perhaps the marginal Pearson correlation coefficient. We follow the vast majority of the neuroscience literature and measure connectivity in each individual by calculating Pearson correlations using samples from a pair of regions. The correlation matrix is then converted into a binary network to represent networks amongst ROIs; see more details in Section 4. In our analysis we have also considered partial correlation matrices (Meinshausen and Bühlmann (2006)) and found that our main results and qualitative findings remain similar.

A number of scientifically important questions are to be addressed for this study. First, which brain regions are activated during these cognitive social tasks and how do these regions function together? Second, if and how subject's social covariates, such as social distress, affect the task-evoked brain connectivity. Third, whether sex differences in brain connectivity during cognitive tasks exist, and if so, how do social covariates influence these differences?

We organize our paper as follows. Section 2 introduces the dynamic network response model and the estimation algorithm. Section 3 presents the simulations, and Section 4 analyzes the task-related study on social cognition and discusses our findings in answering the aforementioned research questions. Section 5 concludes the paper with a short discussion.

Throughout this paper we employ the following notation. Let • denote the outer product and [k] = {1, 2 . . . , k}. For a vector b ∈ R d 1 , let b 2 denote its Euclidean norm. For a matrix B ∈ R d 1 ×d 2 , let B i• and B •j denote its ith row and j th column, B ij denote its (i, j )th entry, and B F denote its Frobenius Norm, respectively. For a tensor B ∈ R d 1 ×d 2 ×d 3 , let B ij k denotes its (i, j, k)th entry, B ij • denote the (i, j )th tube fiber, and B ••k denote the kth frontal slice. For b ∈ R d 3 and B ∈ R d 1 ×d 2 ×d 3 , we define the tensor vector multiplication as

2.1. The dynamic network response model. Consider dynamic networks denoted by G i (V, E i (t)), i ∈ [N], observed from N subjects, where V represents the common set of n nodes and E i (t) represents the set of edges at time point t for subject i. For each subject we also observe a p-vector of covariates, denoted by x i = (x i1 , . . . , x ip ) . At each time point t, the network G i (V, E i (t)) can be uniquely represented by its n × n adjacency matrix A (i) (t), where A (i) jj (t) denotes the edge between nodes j and j at time point t in subject i. The edges can be continuous, binary or nonnegative integers. Without loss of generality, we assume t ∈ [0, 1], and

, where the expectation E(•) is applied elementwise to entries in A (i) (t). We assume that, conditioning on x i , the entries in A (i) (t) are independent and follow an exponential distribution with a canonical link function that

where B 0 (t) ∈ R n×n characterizes the population-level time-varying network connectivity and B l (t) ∈ R n×n characterizes the time-varying effects of the lth covariate on the network connectivity. The function g(•) is an invertible link function, as commonly used in GLMs (McCullagh and Nelder (1989)), and is applied elementwise to entries in μ (i) (t).

Let B ljj (t) denote the (j, j )th element of B l (t). To estimate the unknown functions B ljj (t)'s, we consider a nonparametric estimation using B-spline approximations. Specifically, we approximate B 0jj (t)'s using a K 0 -dimensional basis, denoted by φ 0 (t) = (φ 01 (t), . . . , φ 0K (t)) such that B 0jj (t) ≈ φ 0 (t) × b 0jj , and approximate B ljj (t)'s using a K 1 -dimensional basis, denoted by φ 1 (t) = (φ 11 (t), . . . , φ 1K (t)) , such that B ljj (t) ≈ φ 1 (t) × b ljj . Defining B 0 ∈ R n×n×K 0 and B l ∈ R n×n×K 1 such that B ljj . = b ljj for all j , j and l, model (2) can be rewritten as

where × 3 is defined as in (1), and B 0 ∈ R n×n×K 0 , B 1 . . . , B p ∈ R n×n×K 1 are unknown tensor coefficients. A graphical illustration of model ( 3) is given in Figure 1.

One challenge in estimating model ( 3) is the inherent high-dimensionality of the tensor coefficients. In our analysis of the HCP social cognition study, if we set K 0 = K 1 = 10, then each coefficient tensor B l is of dimension 68 × 68 × 10 = 46,240, far exceeding the number of subjects in the study. Thus, it is imperative to employ effective dimension reduction assumptions that can facilitate estimability and interpretability. Next, we move to discuss the dimension reduction assumptions placed on the baseline effect coefficient tensor B 0 and the covariate effect coefficient tensors B 1 , . . . , B p . We also discuss the need for considering different assumptions for these two types of effects: Low-rankness on B 0 . The component B 0 is the baseline coefficient tensor, and we assume that it possesses a low-rank structure. This specification assumes that there is a lowdimensional structure in the baseline time-varying network connectivity such that both the nodes and the basis coefficients have lower dimensional representations. This is similar to, but more general than, for example, the stochastic blockmodel (Holland, Laskey and Leinhardt (1983)), a well-studied network model that assumes the nodes form a number of groups, and after reorganizing by group membership, the connecting probability matrix is a block matrix.

In our data problem, the low-rank assumption effectively reduces the number of parameters and increases computational efficiency. Specifically, we assume that B 0 admits the following rank-R CP decomposition (Kolda and Bader (2009)):

where w r ∈ R + , u 1r ∈ R n , and u 3r ∈ R K 0 . For identifiability, we assume u 1r 's and u 3r 's are unit length vectors. We note that the above formulation is for undirected networks. When the networks are directed, we can write B 0 = R r=1 w r u 1r • u 2r • u 3r , where u 2r ∈ R n is a unit length vector. Notably, CP low-rank is a special case of the Tucker low-rank (Kolda and Bader (2009)). The definition of the Tucker decomposition for B is expressed as

and U 3 ∈ R r 3 ×d 3 are matrices with orthonormal columns. It is seen that the CP low-rank is a special case of Tucker low-rank with G l 1 l 2 l 3 = 0, unless l 1 = l 2 = l 3 (i.e., superdiagonal). Although the Tucker decomposition is more flexible, it also brings computational challenges, due to the orthonormal identifiability constraint on U 1 , U 2 , U 3 and an increased number of parameters. As a result, the CP low-rank structure is commonly considered in tensor methods for neuroimaging data, such as tensor regressions (Zhou, Li and Zhu (2013), Sun et al. (2017), Zhang and Li (2017), Zhou et al. (2021)), tensor clustering (Cai, Zhang and Sun (2021)) and covariance decomposition (Deng, Tang and Qu (2023)), and has been found to give a good performance.

Structured sparsity in B 1 , . . . , B p . We assume that the subject covariates have sparse effects on the dynamic network connectivity; that is, the effects concentrate on a small number of edges. This is scientifically plausible, as brain connections are energy consuming and biological units tend to minimize energy-consuming activities (Bullmore and Sporns (2009)). Sparsity also greatly reduces the number of free parameters and improves interpretation of the resulting model. Specifically, we assume that B l is structurally sparse in that most tube fibers in B l are zero, corresponding to sparse covariate effects. In particular, tube fiber B ljj • = 0 indicates covariate l has no effect on edge (j, j ), that is, B ljj (t) = 0 for all t. To encourage structural sparsity, we consider the group lasso (Yuan and Lin (2006)) penalty, defined as ( 4)

Different assumptions on B 0 and B 1 , . . . , B p . We briefly discuss the benefits and necessity of imposing separate structures on B 0 and B 1 , . . . , B p . It is natural to think that one could stack B 0 , B 1 , . . . , B p into one higher-order coefficient tensor of size n × n × K × (p + 1) (assume K 0 = K 1 = K) and specify it to be both low-rank and sparse. However, assuming B 0 to be sparse may not be plausible in the GLM setting. For instance, when the network edges are binary and g(•) is the logit link, g(0) yields a connecting probability of 0.5; when the network edges are counts and g(•) is the log link, g(0) is not well defined. Correspondingly, a sparse B 0 does not necessarily imply sparsity in the baseline connectivity and may not even be well defined. This issue is unique in using sparse GLM to model edges in a network. Finally, more complex structures on B 1 , . . . , B p can be incorporated (e.g., B 1 , . . . , B p are low-rank and sparse), which can further reduce the number of effective parameters. However, such assumptions are expected to incur a much higher computational cost and also involve more tuning parameters on, for example, the rank of each coefficient. To balance model complexity and feasibility, we focus on the current assumption that assumes B 1 , . . . , B p have structured sparsity. 3) the negative log-likelihood function, up to a constant, can be written as

where McCullagh and Nelder (1989)). We estimate the parameters w, U 1 , U 3 , and by solving the following optimization problem:

where P(•) is as defined in (4) and λ is a tuning parameter. The optimization problem in ( 6) is computationally challenging, as the size of the networks, the dimension of the covariates, and the number of basis functions can be large in practice. The GLM loss function further increases the computation burden due to its nonlinearity. While (6) is nonconvex, the conditional optimization with respect to u 1r , while fixing all other parameters, is convex, and the same holds for w, u r3 's, and B j 's. This observation permits an alternating minimization algorithm. One potential issue in such an approach is that solving for , conditional on all other parameters, is a regularized optimization problem of dimension n × n × K 1 × p. This can be computationally expensive when the network size n, the number of splines K 0 , K 1 , and the dimension of the covariates p are large. To tackle this challenge, we consider a proximal gradient descent algorithm that is easy to implement and computationally efficient. Our estimation procedure is summarized in Algorithm 1.

In Step 2, ũjr 's are solved using a Newton-type algorithm (Schnabel, Koonatz and Weiss (1985)), and the gradients are given in Section S7.1 in the Supplementary Material. In Step 3

Algorithm 1 Optimization procedure of ( 6)

Input: rank R, tuning parameter λ, bases dimensions {K 0 ,K 1 } and step size η.

Step 1: initialize w (0) , U (0)

p . Repeat Steps 2-5 for t = 0, 1, . . . until convergence.

Step 2: repeat the following steps for r = 1, 2, . . . R. ũ(t+1)

3R , (t) ).

Step 3:

3

),

Step 5: repeat the following steps for s = 0, 1, . . . until convergence.

(t,s+1) = S λη ( (t,s) -η∇ ℓ(w (t+1) , U

we define two matrix operators for U . Norm(U ) calculates the ℓ 2 norms of columns in a matrix U , and Unit(U ) rescales the columns of a matrix into unit vectors. That is,

In Step 5 we employ the fast iterative shrinkage-thresholding method (FISTA, Beck and Teboulle, 2009) under group lasso penalty. Specifically, we define the shrinkage operator by S λη ( ) = (T λη (B 1 ), . . . , T λη (B p )) ∈ R n×n×K 1 ×p , where

and (x) + = max(0, x). In the FISTA algorithm and at step s + 1, the iterative shrinkage operator S λη (•) is not directly applied to the previous point (t,s) but rather at the point (t,s) , which uses a specific linear combination of the previous two points (t,s) and (t,s-1) . The FISTA algorithm has been shown to enjoy a fast global rate of convergence (Beck and Teboulle (2009)) and is easy to implement. The stepsize η is typically chosen as the Lipschitz constant of ∇ ℓ(w, U 1 , U 3 , ), which can be approximately calculated given the initial values.

In Algorithm 1 we need to determine the initial values for the alternating minimization procedure. To obtain a good initial estimate, we first estimate B (0) 0 , B (0) 1 , . . . , B (0) p via an elementwise generalized spline regression; see (8). We then estimate w (0) , U

1 , U (0) 3 via a CP decomposition of the estimated B (0) 0 . In our experiments this initialization procedure leads to a good numerical performance of Algorithm 1. The accuracy of this initialization procedure is evaluated in Section 3.

Parameter tuning. The B-spline dimensions K 0 and K 1 , rank R, and regularization parameter λ are tuning parameters in our algorithm. We choose these parameters using the eBIC criterion that was first developed for variable selection in the diverging dimension regime in Chen and Chen (2012). It has been demonstrated that the eBIC function is effective as a heuristic criterion to balance model fitting and complexity when used in low-rank estimation problems (Srivastava, Engelhardt and Dunson (2017), Cai, Zhang and Sun (2021), Zhang, Sun and Li (2023)). Specifically, we choose the combination of (K 0 , K 1 , R, λ) that minimizes

where ℓ is the loss function in (5), and ŵ, Û 1 , Û 3 , ˆ are the estimates of w, U 1 , U 3 , under the working rank and regularization parameter. In our numerical experiments, the above eBIC is found to be minimized at the true rank and sparsity level under the selected λ (see Section S3). To reduce computational cost in tuning, we recommend using cubic B-splines in practice. If desired, the order of the B-splines can be another tuning parameter selected using the eBIC function.

We conduct simulations to investigate the performance of our proposed method. We focus on symmetric networks and compare our proposed dynamic network response regression method, referred as DNetReg, with two alternative elementwise approaches. To reduce computational cost, we set K 0 = K 1 = K, and the spline dimension K is set to the true value in simulations; see Section S3 of more details of tuning K, R, and s 0 using the eBIC.

The first elementwise approach, referred as EdgeReg, fits elementwise GLMs at each time point t k . That is, for any j, j ∈

This elementwise approach ignores both the network structure and the temporal smoothness in the dynamic brain connectivity. The second elementwise approach, referred as DEdgeReg, fits a generalized spline regression to each entry in A jj (t). Specifically, for any j, j

A Newton-type algorithm is employed to estimate the parameters in the above model. The method DEdgeReg is used to find the initial values in Algorithm 1. We simulate N binary dynamic networks of size n × n in [0, 1] from model (3), where A jj (t), t ∈ [0, 1], follows a Bernoulli distribution and g(•) is taken to be the logit link function. The covariates x i 's are generated independently from N (0, 1), and we standardize the columns of the design matrix to have zero mean and unit standard deviation. For B 0 = R r=1 w r u 1r • u 1r • u 3r , we first generate the entries of u 1r and u 3r from N (0, 1), set w r = u 1r 2 2 u 3r , and then we standardize u 1r and u 3r as unit length vectors. For B 1 we randomly set s 0 proportion of its entries to be 1 and the rest to zero such that s 0 = B 1 0 /(n 2 K 1 ). The basis functions in φ 0 (t) and φ 1 (t) are set to B-spline bases with K 0 = K 1 = 8, and the knots are equally spaced in [0, 1].

To evaluate the estimation accuracy, we report estimation errors B 0 -B0 F , B 1 -B1 F , and N i=1 μ (i) -μ(i) F /N , where μ(i) (t) = g -1 ( B0 × 3 φ 0 (t) + x i ( B1 × 3 φ 1 (t))). Furthermore, to evaluate the edge selection accuracy from our method, we report the true positive rate (TPR) and false positive rate (FPR) in identifying the nonzero entries in B 1 . The first elementwise approach EdgeReg does not estimate spline coefficients B 0 and B 1 , and thus their estimation errors are not reported. While estimates from EdgeReg are not sparse, the p-values for B ljj (t h )'s are directly available from standard GLM model fitting. In our evaluations we apply Bonferroni correction to these p-values and then calculate the TPR and FPR in identifying the edges modulated by x 1 , that is, entries (j, j )'s with nonzero time-varying covariate effects B 1jj (t)'s. Specifically, we define P BC ∈ R n×n×T , where P BC jj h is the pvalue in evaluating the significance of B 1jj (t h ) from ( 7), after the Bonferroni correction of n × n × T tests. Defining H ∈ R n×n with H jj = 1{min(P BC jj . ) ≤ 0.05}, and

The FPR and TPR are calculated as

where * denotes the elementwise product. The second elementwise approach DEdgeReg does not give sparse estimates, and there are no readily available inference results to calculate pvalues; hence, their TPRs and FPRs are not reported.

We set the number of subjects N = 50, 100, the number of equally spaced time points T = 100, and consider the number of nodes n = 50, 100, rank R = 2, 5, and the sparsity proportion s 0 = 0.05, 0.1, respectively. Tables 1 and 2 report the average accuracy measures over 50 replications with sample size N = 50, 100, respectively, with the standard deviations shown in parentheses. It is seen that our proposed method achieves the best performance among all competing methods, in terms of both estimation accuracy and selection accuracy, and this holds for different sample sizes N , numbers of nodes n, ranks R, and sparsity levels s 0 . Moreover, the estimation error of our method DNetReg decreases as network size n, rank R, and sparsity proportion s 0 decrease and as sample size N increases. Estimation errors from EdgeReg and DEdgeReg are not sensitive to R or s 0 , as they are elementwise approaches and do not consider the low-rank and sparsity structure in the tensor coefficients. In terms of edge selection accuracy, EdgeReg is overly conservative after the Bonferroni correction, and its TPRs are close to zero. In our analysis we also considered FDR (or BH) correction (Benjamini and Hochberg (1995)) for p-value corrections, and the results are similar. We further report the eBIC over varying R, K 0 , K 1 and s 0 values in Section S3 of the Supplementary Material, consider sensitivity analysis of our method under model misspecifications in Section S4 of the Supplementary Material, and discuss the computational cost of DNetReg in Section S5 of the Supplementary Material.

We consider the social cognition study in the HCP study, a publicly available HCP study with preprocessed cortical-surface fMRI data (see https://www.humanconnectome.org/study/hcp-young-adult/document/900-subjects-data-release). 2 In our analysis we mapped the preprocessed data to the Desikan-Killiany Atlas (Desikan et al. (2006)). The social cognition study collects task-related fMRI data from N = 843 healthy adult subjects. Specifically, the fMRI data are collected on 274 evenly spaced time points covering an initiation countdown (five seconds) followed by five video blocks (23 seconds each) with fixation blocks in between (15 seconds each). The first 11 scans in the initiation countdown period are removed in our analysis. The fMRI data are then preprocessed and summarized as a 68×263 spatial-temporal matrix for each subject using the Desikan-Killiany Atlas (Desikan et al. (2006)) with n = 68 ROIs (see Table S4). For each subject the dynamic network is constructed by calculating a sequence of connectivity matrices of dimension 68 × 68 over T sliding windows, each summarizing the connectivity between the 68 brain regions in a given window. We let the number of samples in each window and the overlap between adjacent windows be 30 and 5, respectively, giving a total of T = 47 networks per subject. To check whether the results are sensitive to the parameters used to construct the dynamic connectivity matrices, we also consider setting the window size to 25 and 35 and step size 3 and 7. The result are shown in Section S6.1 of the Supplementary Material. We determine connectivity in each individual by computing Pearson correlations between samples from a pair of regions and create binary networks by setting A jj (t h ) = 1 if the computed correlation value is greater than 0.5 and A jj (t h ) = 0 otherwise, and this gives an average network density about 15%. This procedure can eliminate weak functional connectivity and is commonly employed in existing neuroscience literature (Power et al. (2011)). In our analysis we have also considered partial correlation matrices (Meinshausen and Bühlmann (2006)) and applied other thresholding values, such as 0.6, to the Pearson correlation matrix and found that our main results and qualitative findings remain similar. Notably, we have also considered continuous edge weights, calcu- 3) with Gaussian distribution and g(•) as the identity link, and the results are shown in Section S6.2 of the Supplementary Material. Our analysis is carried out in R, and the code is available online. 3In the social cognition study, there are 374 males and 469 females, aged between 22 and 36 years old. In addition, the HCP data contain social covariates on three categories, including friendship and loneliness in the Companionship category, perceived hostility and perceived rejection in the Social Distress category, and emotional support and instrumental support in the Social Support category. Our preliminary analysis finds that these covariates are highly correlated with correlations ranging from 0.4 to 0.6 and a median of 0.5. Due to this consideration and to reduce computational cost, we focus our analysis on perceived hostility (e.g., how often people argue with me, yell at me, or criticize me) representing the Social Distress category. A higher perceived hostility shows increased social distress, which is the extent to which an individual perceives his/her daily social interactions as negative or distressing (Lieberman (2007)). Moreover, we have also considered the friendship covariate from the Companionship category and the emotional support from the Social Support category. The model fitting results are shown in Section S6.3 of the Supplementary Material. The goal of our analysis is to characterize the baseline brain connectivity in tasks, to ascertain how social covariates modulate the subject-level connectivity changes, and to examine whether there are any sex-specific differences. We apply our proposed model to the dynamic connectivity networks from males and females, respectively (see more discussions in Section S6.3). The social covariate is standardized to have mean zero and variance one, and we consider cubic B-spline basis with different equally spaced knots of B 0 and B 1 . Using the eBIC function, the basis dimensions were selected as K 0 = 8, K 1 = 10, the rank as R = 9, and the sparsity proportion as s 0 = 0.11 for males and K 0 = 12, K 1 = 10, R = 11, and s 0 = 0.17 for females. The results of fitting the model with cubic natural spline are shown in Section S6.4 of the Supplementary Material.

Baseline brain connectivity. We start by examining the estimated baseline connectivity coefficient B0 . Figure 2 plots the baseline connectivity averaged over time, that is, g -1 ( T h=1 B0 × 3 φ 0 (t h )), where g(•) is the logit link function and nodes are organized by results from a K-means clustering. Specifically, we apply K-means clustering, based on SVD of the average connectivity matrix T h=1 B0 × 3 φ 0 (t h ) for male, and identify five clusters among the 68 ROIs. The members of each cluster are given in Table 3. While clustering results using B 0 estimated for females are similar, we use the same clustering labels to facilitate comparisons. Anatomically, the first community contains mostly nodes in the cingulate gyrus, the second and fifth communities contain nodes from the frontal lobe, the third community contains nodes from the temporal lobe, and the fourth community contains nodes from the frontal, parietal, occipital, and temporal lobes (see Tables 3 and S4). Many of the 68 anatomic ROIs in the Desikan Atlas overlap with the resting-state functional modules. We find that community 1 is associated with emotion formation and processing, commu-

TABLE 3 The anatomic regions of interest in the identified communities 1 Caudalanteriorcingulate, isthmuscingulate, paracentral, posteriorcingulate, transversetemporal, insula 2 Cuneus, lingual, pericalcarine, postcentral, precentral, precuneus, rostralmiddlefrontal, superiorfrontal, supramarginal 3 Entorhinal, parahippocampal, temporalpole 4 Bankssts, caudalmiddlefrontal, fusiform, inferiorparietal, inferiortemporal, lateraloccipital, middletemporal, parsopercularis, parstriangularis, superiorparietal, superiortemporal 5 Lateralorbitofrontal, medialorbitofrontal, parsorbitalis, rostralanteriorcingulate, frontalpole nity 2 is related to visual, attention, and emotion regulation modules, and community 4 is enriched with visual and object identification. The lateral occipital gyrus in community 4, lingual gyrus in community 2, and pericalcarine gyrus in community 2 are from the occipital lobe, a region responsible for interpreting the visual world (Goldenberg et al. (1991)), and is seen to be active for both males and females. For both males and females, we find that connectivity between communities 2 and 4 is more active both within and between the two hemispheres, especially the temporal parietal junction, superior temporal cortex regions, and occipital gyrus, which are all relevant in social cognition. This is in line with previous research which showed that mental animations stimulate these regions (Castelli et al. (2000), Barch et al. (2013)). Within each hemisphere, males have higher connectivity within communities 2 and 4, and this is consistent with the existing findings that males have increased intrahemispheric connectivity (Ingalhalikar et al. (2014)).

Social effects on brain connectivity and sex differences. We next examine the estimated covariate effect coefficient B1 . Figure 3 plots the heatmap of estimates for males and females, where we show ( T t=1 B1 × 3 φ 1 (t))/T , the mean coefficient of B1 , representing the covariate effect on brain connectivity averaged over time.

It is seen that the social effects on connectivity show different patterns in males and females. Specifically, the estimated B1 has sparsity portions equal to 0.17 and 0.11 for females and males, respectively. Hence, the social effect on connectivity is more sparse in males, and such differences are observed in within-and between-community connectivity within and across hemispheres. Compared to males, the social covariate is seen to more notably decrease the connectivity between communities 2 and 4 within the right hemisphere and also across hemispheres in females, suggesting that the task-related brain connectivity in females is more sensitive to social stress. This supports existing findings that social stress influences brain connectivity and emotional perception differently for males and females (Mather et al. (2010)). In general, the perceived hostile social distress covariate has a negative impact on the connection response for females both within and between communities, particularly for community 4, while it tends to have a positive impact on the connection response for males. The above findings on sex-specific difference are interesting, and they may be linked to existing research on sex differences in neural response to psychological stress (Wang et al. (2007)).

Figure 4 shows the social effects on brain connectivity in males and females during different periods of the experiments including watching a mental video, resting, and watching a random video. It is seen that, during a mental video, the connectivity between temporal and occipital lobes in females is more affected by social stress. The temporal lobe plays an important role in visual perception and processing emotions, and the occipital lobe is related to visual processing, containing most of the anatomical region of the visual cortex (Goldenberg et al. (1991)). Moreover, We discover that females have greater activity in the across-lobe connectivity, particularly among the temporal, parietal, and occipital lobes (Ingalhalikar et al. ( 2014)) in all periods. This finding suggests some interesting patterns that warrant further investigation and validation.

Finally, we have visualized elementwise plots of estimated baseline effect B 0 (t) and covariate effect B 1 (t) in Figure 5 below. Our presentation focuses on region 22 in baseline effect B 0 (t) and regions 10 and 28 in covariate effect B 1 (t), as they show interesting differences between females and males. Specifically, plots (a)-(b) show the time-varying baseline effects in B 0 (t) between the left posterior cingulate (region 22) and other regions. For females, there are five notable valleys near time points t = 5, 11, 20, 32, and 42, corresponding precisely to three mental tasks and two random video clips. In contrast, these trends are much less noticeable in males, suggesting the posterior cingulate is more responsive in females during video viewing. This is supported by existing findings that the posterior cingulate gyrus is more active in females during emotion-related tasks, given its role in emotion regulation and processing (Proverbio et al. (2009)). Plots and females, the connectivity between Region 62 and left superior parietal (28) is the most affected by the social covariate, reflecting the important role of the superior parietal lobule in cognitive tasks (Alahmadi (2021)). It is seen that, for females, the connectivity between Region 62 and others are more affected by the social covariate.

A permutation based procedure to examine sex differences. Developing the asymptotic distribution of the estimated B 1 under the CP low-rank and sparsity constraints in our model is challenging. In this section we conduct an ad hoc permutation based procedure to examine whether the previously identified sex-specific differences are meaningful.

Specifically, we randomly permute the sex labels across subjects 100 times. In each permutation i, we divide the N = 843 samples into two groups, based on the permuted sex labels, and apply the proposed model to the male and female groups, respectively. We denote the coefficient tensors as B male,i 0 (or B female,i 0 ) and B male,i 1 (or B female,i 1 ) in permutation i, i ∈ [100]. To quantify the difference in B 1 between males and females, we calculate the absolute distance between the coefficient vectors for each (j, j ). Specifically, we write We define a binary matrix S ∈ R n×n S jj = 1 100 i=1 1 D obs jj > D per,i jj ≥ 95 , where 1(•) is the indicator function. Correspondingly, S jj = 1 if the observed sex difference is the same as or greater than the 95th percentile of permuted sex difference. Figure 6(b) plots S, which further illustrates that the sex differences within community 4 and between communities 2 and 4 are likely significant (regions in the blue and black dashed lines), affirming the findings in Figure 3. We also consider comparing results based on subgraphs of interests, shown in Figure S5, where sex-specific differences from observed data are consistently greater than those from permuted data. FIG. 7. Heatmaps of ( T t=1 B1 × 3 φ 1 (t))/T estimated by DEdgeReg, with rows and columns ordered the same as Figure 4. The top and bottom panels show the results without and with thresholding, respectively. 4.2. Results using existing methods. We evaluate the performance of two alternative methods including an elementwise method DEdgeReg, evaluated in Section 3, and GLSNet (Zhang, Sun and Li (2023)), a nontime-varying matrix response regression model. Since GLSNet is not designed to model dynamic networks, we directly calculate the connectivity matrix based on all 263 scans using the same procedure that binarizes the Pearson correlation matrix. Using GLSNet and the recommended eBIC function in Zhang, Sun and Li (2023), the rank is selected as R = 5 and the sparsity proportion as s 0 = 0.025 for males, and R = 13 and s 0 = 0.0359 for females.

Figure 7 shows ( T t=1 B1 × 3 φ 1 (t))/T (representing the effect averaged over time) estimated by DEdgeReg with or without threshlding at ±0.1. It is seen that the estimates from the elementwise method DEdgeReg are very noisy, and they identify a large number of regions with relatively small signals. The estimated social score effect coefficients B1 from GLSNet are shown in Figure 8. For both males and females, the estimates are highly sparse. In males, several areas associated with social cognition, such as the temporal parietal junction, superior temporal cortical regions, and occipital gyrus, do not appear to be engaged. This can potentially due to the fact that GLNet ignores the dynamic changes of brain connectivity during the experiments.

In this paper we study the task-evoked brain connectivity by introducing a new semiparametric dynamic network response regression that relates a dynamic brain connectivity network to a vector of subject-level covariates. A key advantage of our method is to exploit the structure of dynamic imaging coefficients in the form of high-order tensors. We briefly comment on potential future research. In our model setup, we assume that the tensor coefficients B 1 , . . . , B p are sparse. More complex structures, such as the low-rank or fused structures, can be considered as well, though they will increase the computation time and complexity in tuning. We have included two extended models with time-varying covariates and low-rank covariate effects. The corresponding simulation results and real data analyses are presented in Sections S1 and S2 of the Supplementary Material, respectively. In Section 4.1 we consider an ad hoc permutation procedure to evaluate the identified sexspecific differences. A more rigorous approach would be to derive the asymptotic distribution of B 1 and carry out hypothesis testing. This is not a trivial task due to the involvement of both low-rank and sparse constraints on the model parameters. We leave this investigation to future research.