Role of the EEG Theta Network During Software Production: A Connectivity Study

Software programming is an acquired evolutionary skill originating from consolidated cognitive functions (i.e., attentive, logical, coordination, mathematic calculation, and language comprehension), but the underlying neurophysiological processes are still not completely known. In the present study, we investigated and compared the brain activities supporting realistic programming, text and code reading tasks, analyzing Electroencephalographic (EEG) signals acquired from 11 experienced programmers. Multichannel spectral analysis and a phase-based effective connectivity study were carried out. Our results highlighted that both realistic programming and reading tasks are supported by modulations of the Theta fronto-parietal network, in which parietal areas behave as sources of information, while frontal areas behave as receivers. Nevertheless, during realistic programming, both an increase in Theta power and changes in network topology emerged, suggesting a task-related adaptation of the supporting network system. This reorganization mainly regarded the parietal area, which assumes a prominent role, increasing its hub functioning and its connectivity in the network in terms of centrality and degree.


I. INTRODUCTION
E VERYDAY life and work activities require the brain to perform complex functions, such as paying attention, acquire and process complex information, design targeted responses, plan and execute actions. Neuroergonomics is a branch of the cognitive neuroscience research that focuses on the study of the human brain in everyday-life scenarios. One of its goals is the implementation of new technologies to optimize the efficiency and productivity of workers by continuously monitoring their physiological responses and mental state [1], [2]. A popular application lies in the context of software engineering, in which, nowadays, industries worldwide are investing [3], [4], [5].
The idea is that by understanding the neural mechanisms supporting software development processes, it could be possible to recognize error-prone scenarios and, consequently, to provide code optimization strategies and minimize the cost of software faults (bugs) [5]. It is worth noting that software development represents one of the largest industrial sectors in the world, not only because software development companies are among the biggest companies, but also due to the transversal nature of software, which makes software a necessary element in many other activity sectors. Even a modest improvement in the efficiency and quality of software production will have a huge impact in society. Nevertheless, poor and scattered literature exist about this specific topic [6]. Software programming is a cognitive skill recently learned by humans and its neural basis are still not completely explored and understood. Most of the existing studies are based on functional magnetic resonance (fMRI) and investigate the brain areas involved in software development. Usually, a single activity is studied, such as code inspection and comprehension [6]. A common finding is the involvement of fronto-parietal areas during programming tasks, but some differences in activation patterns are reported depending on the specific task and experiment performed [7], [8], [9], [10]. Siegmund and colleagues showed, for example, that code comprehension tasks, contrasted with syntax tasks, recruit the same fronto-parietal and left lateralized areas associated to working memory and language system, respectively [11], [12]. Krueger and colleagues [7] investigated the differences in brain mechanisms supporting prose and code writing, concluding that this latter is associated to a right lateralization, while prose production elicits a prominent activation of left brain areas associated to the language system. Accordingly, Ivanova et al. [8] suggested a distinction between networks supporting software programming and the ones sustaining language processing. Specifically, the authors observed that code comprehension does not involve the language system but elicits a bilateral activation of fronto-parietal domaingeneral network, supporting working memory and cognitive control. Hence, the most acknowledged hypothesis is that programming abilities rely on pre-existing networks, historically activated for more ancient cognitive tasks, and specialized over the years. Software programming, indeed, involves different functions, such as working memory, problem-solving, language (natural or code) comprehension and production, mathematic calculation, and logic reasoning. Therefore, the recruitment of the same or partially overlapping networks supporting all these cognitive skills could be expected [9]. Nevertheless, trying to describe the brain mechanisms underlying the software programming process considering distinct assignments, as done in the recent literature [6], could lead to an oversimplification. Indeed, software programming in real settings includes different tasks i.e., instructions reading and comprehension, problem solving strategies, code writing and debugging. It is therefore reasonable to expect a dynamic recruitment of different brain networks [8].
To verify this hypothesis, we propose the use of Electroencephalography (EEG) to monitor programmers' brain activity during a realistic programming assignment, while programmers are effectively developing code in a usual software development setup (i.e., not just reading/ comprehending code as in most studies in the literature [6]).
The majority of the few existing works investigating EEG signals during code writing tasks are based on the extraction of standard pre-existing features for the classification of outcome measures, such as task difficulty, programmers' proficiency and cognitive workload [13], [14], [15], [16]. The main EEG-based findings revealed that both Theta rhythms and fronto-parietal networks are involved during programming tasks, suggesting the recruitment of the working memory (WM) system [17], [18].
Therefore, for the first time, to the best of our knowledge, we present the exploration of brain connectivity reorganization in the Theta band in function of different cognitive tasks: natural language text reading, code reading and programming (i.e., code writing). In order to obtain an experimental experience as similar as possible to real-life programming scenarios, in the present study EEG traces were acquired during an event-free realistic programming task.
In literature, a plethora of methodologies for the study of EEG based brain connectivity exist, which are based on different models and mathematical assumptions [19]. A first distinction is between functional, which measures the strength of the reciprocal connection among brain areas, and effective connectivity, which also adds an information of directionality. Effective connectivity can be explored by means of either a multivariate or pairwise approach. The former can better reconstruct complex connectivity pattern, but it implies an increased model complexity [20] and difficult interpretation. Therefore, when a large number of network nodes are recruited, the latter approach can be preferred. Even focusing on a pairwise effective connectivity approach, the literature offers a large variety of methods to choose from [21]. For example, some information theory based metrics, such as the Transfer Entropy [22], offer the possibility to explore nonlinear brain directional couplings. Concerning linear methods, the Granger causality (GC) is an example of widely used method to explore pairwise effective connectivity [19], despite its vulnerability to volume conduction (VC) artifacts. For the present study, the Phase Slope Index (PSI) [23] was chosen among all the available possibilities. It is a pairwise effective connectivity measure that has been proven robust against VC. PSI, indeed, is computed starting from the imagery part of coherence, with the aim of isolating only time-lagged interactions between two time series. Therefore, the instantaneous synchronizations mainly due to VC are neglected.
The present paper is divided in the following sections. Section II, Materials and Methods, describes the considered population, the experimental protocol adopted, and the methods used for the spectral and connectivity analysis. Section IV, Results, summarizes the outcomes of the statistical analysis performed on spectral and connectivity metrics. Section V, Discussion, provides an interpretation of obtained results. Finally, conclusions, remarks and future developments are presented in Section V.

A. Subjects
Electroencephalographic traces were collected from 11 healthy subjects with experience in software programming. Participants were selected among researchers and students at Politecnico di Milano and Università Statale di Milano. Their eligibility and proficiency level were determined by means of an ad-hoc designed screening questionnaire. This latter consisted of multiple-choice questions about the functioning of ten code snippets written in C language. To be part of the study, volunteers had to correctly solve at least four questions out of ten. Specifically, three females and eight males, with mean age 27.72±6.06, were enrolled, 10 of them of native Italian language, 1 of native Chinese language, but all of them were fluent in English. All subjects signed an informed consent, and their data were anonymized. The study was approved by the Ethics Committee of Politecnico di Milano.
As a reward for participation, all subjects received a voucher of 40 euros for online purchases at the end of the experiment.

B. Experimental Setup and Protocol
Physiological signal acquisitions were performed with the SD LTM 64 express polygraph recording system (Micromed, Mogliano Veneto, Italy), which uses an elastic EEG cap with 61 electrodes placed according an extension of the 10-20 international system [24]. EEG traces were acquired with a sampling frequency of 256 Hz in monopolar configuration, with the reference electrode placed between CPz and Pz.
During the whole experiment subjects were sit in front of a PC screen, in a relaxed, dimly lit and silent environment [17].
The experimental protocol was organized in two runs, and in each of them, participants performed three tasks: (1) natural text reading, (2) code reading and (3) realistic programming. These three tasks were presented in random order and interspersed with a resting phase of 30 seconds, during which a screen with a cross in the middle was presented to the volunteer. The scheme of the described protocol is shown in Fig. 1.
Specifically, during the natural text reading task, volunteers had to read a text written in English language for 60 seconds. The code reading task lasted 5 minutes and consisted in the reading and comprehension of a code snippet written in C language.
During the programming phase, volunteers had to develop a code in C language following the instructions provided at the beginning of the assignment. This task had a minimum duration of 5 minutes, after which volunteers were allowed to end the exercise any time before the maximum duration of 20 minutes.
At the end of each run, volunteers compiled an evaluation form in which they were asked to rate the programming task in terms of mental effort, task fulfillment, pressure over time and discomfort.
The experiment was managed by means of a code in-house written in MATLAB, while volunteers used the Eclipse IDE for the programming tasks.

C. EEG Pre-Processing
The whole pre-processing pipeline was performed by means of the EEGLab toolbox [25] and in-house MATLAB scripts.
Raw EEG signals were filtered with a zero-phase band-pass FIR filter in the frequency range 1-45 Hz and then downsampled at 128 Hz. The order of the FIR filter was set to ensure a suppression of -6 dB at the cut-off frequencies. After a visual inspection, noisy segments and bad channels were removed from EEG traces. Then, Independent Component Analysis (ICA) through Infomax algorithm [26] was executed. The EEGLab plug-in ICLABEL [27] was used to support the identification of artifactual sources [28]. Specifically, this plug-in attributes to each component the probability of belonging to seven main source classes: i) Brain, ii) Muscle, iii) Eye, iv) Heart, v) Line Noise, vi) Channel Noise, and vii) Other.
After the removal of artifactual components, bad channels interpolation was performed.
As a final step, a modified Common Average Reference (CAR) was applied: the average signal was computed on the 19 electrodes of the International 10/20 System and then removed from all the 29 acquired traces to avoid possible unbalancing due to the frontal lobe spatial oversampling.

D. Spectral Analysis
Thirty-second-long artifact free EEG segments were extracted from the fixation phase (FIX), the natural text reading and (NTR), code reading (CR) tasks, while three different segments of equal length were extracted from the programming task. Precisely, EEG segments at the start (P1), in the middle (P2) and at the end (P3) of the task were considered. The stationarity of the considered segments was tested performing the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test, using a p-value of 0.01 [29], [30].
The power spectral density of each segment was estimated using the Welch's method, with 1-second-length windows hamming tapering and 50% overlap. EEG signals in each window were normalized to zero mean and unitary variance. Then, the relative power in the standard Theta frequency bands (4-8 Hz) was computed for each electrode. Finally, the percentage power variation with respect to the resting condition (i.e., FIX) was computed for each task according to equation (1): where, P task is the power in the Theta band during the considered segment (i.e., NTR, CR, P1, P2 or P3), whileP F I X is the average power in the Theta band during FIX. Nine brain areas were defined as clusters of EEG channels (TABLE I)

E. Connectivity Analysis
Pairwise effective connectivity (EC) relations between each couple of EEG channels were estimated using the Phase Slope Index (PSI) [23].
The PSI is derived from the imaginary part of the complex coherence and estimates the strength and the direction of information flux between two time series i and j as in equation (2) where C i j is the complex coherence between electrodes i and j, i as reference, C * i j its complex conjugate, ℑ (·) is the operator that isolates the imagery part, δ f is the frequency resolution, and F is the frequency band considered.
The coherence between pairs of EEG signals was computed according to equation (3).
where S i j is the cross-spectrum between i and j, while S ii and S j j are their auto-spectra. These were computed according to the procedure described in paragraph II-D, i.e., by the application of the Welch's method using windows of 1 s with 50% overlap, leading to a frequency resolution of 1 Hz. Such window length allows to properly study brain mechanisms in the Theta band, since assures the detection of at least 4 cycles of the considered frequencies (a fundamental frequency of 4 Hz is, indeed, associated to a period of 0.25 s). Moreover, it also allows to meet an acceptable tradeoff among frequency resolution, variance associated to spectral estimation, robustness of phase difference estimation and the number of values averaged to compute the final PSI in the frequency range F. The obtained connectivity matrix is antisymmetric (i.e.,ψ i j = −ψ i j ): positive values of PSI indicate an information flow that goes from i to j (i.e., rows to columns), while negative values are associated to the opposite direction. Moreover, the net flux index was computed for each node starting from the PSI matrix as in equation (4): This metric considers all the inward and outward connections of each node and indicates, on average, if it behaves as a source (i.e., positive net flux) or as a sink (i.e., negative net flux) of the network.

F. Network and Graph Analysis
As a first step, transformed connectivity matrices were obtained setting at zero non-significant (Fig. 2B) and negative values (Fig. 2C) and applying a final binarization (Fig. 2D) Specifically, negative PSI values were set at zero in order to remove redundant information. In this way, an antisymmetric matrix of positive values is finally obtained, where the directional information is maintained (i.e., information flux directed from rows to columns). Significant connections (PSI values between pairs of signals) were assessed performing a surrogation test. For each couple of channels, 100 couples of surrogates were obtained through the well-known phase randomization method [31] and, for each of them, the PSI was computed. Specifically, the Fourier Transform (FT) method was used: i) the power spectra of each time series were computed by means of the FT, ii) the modulus of the power spectra was maintained, while the phases were randomized with values in the range [−π ; +π ] 100 times, iii) the 100 surrogates were obtained by means of the inverse FT. For each pair of channels, a null distribution of 100 PSI values was generated by computing the index at each iteration of the surrogation procedure. Then, the 95th percentile of such distribution was used as significance threshold: all the connectivity measures that were, in absolute value, lower than this threshold, were set to zero. In the end, matrices of significant PSI values were binarized, that is, each significant ψ i j (connection directed from i to j) were set to 1, the non-significant ones were set to 0. Since the final matrix is positive and asymmetric, the information of directionality is still maintained, i.e., a value equal to 1 in correspondence of row i and column j identifies a significant connection that goes from i to j. Finally, the obtained transformed matrices were used for the computation of reduced networks (RN) and the performance of Graph theoretical analysis. 1) Reduced Networks: The transformed matrices were further manipulated to obtain reduced networks (RNs): starting from the 29 × 29 binarized matrix, 9 brain regions of interest were identified and the significant connections among them were further considered. Connections within each region were discarded. Precisely, the nine regions are the ones listed in TABLE I. The final result is a 9 × 9 connectivity matrix, in which each node of the network is associated to a brain region, while the links among them are obtained by averaging the significant (1) and non-significant (0) connection values between the nodes belonging to each pair of regions. In this way, the connectivity links were normalized between 0 and 1. Moreover, in the computation of such matrices, the links between couples of frontal regions and between the parietal and the occipital regions were neglected because of their spatial proximity. Reduced networks of run 1 and run 2 were averaged for each subject.
2) Graph Analysis: The graph analysis was carried out on the binarized 29 × 29 connectivity matrices. Specifically, each connectivity matrix was associated to a graph, where the nodes and the links were represented by the electrodes and the connectivity measures, respectively. Then, different metrics describing the network topology were computed by means of the Brain Connectivity Toolbox (BCT) [32].
Betweenness centrality (BC), Clustering Coefficient (CC), In Degree (ID) and Out Degree (OD) are the graph measures computed for each node and connectivity matrix. Moreover, the Global Efficiency index (GE) was computed as a whole network measure of integration [33]. BC measures the tendency of each node to behave as a hub of the network and is computed according to the following formula (5): where n is the number of nodes of the network, p l, j (i) represents the shortest path between node l and node j passing through node i and p l, j is the number of shortest paths between nodes i and j. CC indicates of how much the vertices of the graph are organized in functional subnetworks and it is computed as in equation (6).
where t i is the number of triangles around node i, and k i is the degree of the node. i.e., the number of links to which i is connected to. ID and OD count the amount of inward and outward connections for each node and are defined as in (7) and (8), respectively: where a ji and a i j are the binary values associated to the link that goes from node j to node i and to the link that goes from node i to j, respectively. Finally, GE was computed according to (9): where d i j is the shortest path between node i and node j. Graph measures obtained for run 1 and run 2 were averaged. Consequently, average graph measures for FIX, NTR, CR, P1, P2 and P3 were obtained for each subject.
Node-based graph measures were then averaged on the brain areas listed in TABLE I.

G. Statistical Analysis
Non-parametric statistical tests were performed given the small sample size (i.e., 11 subjects).
The statistical significance of the Theta power variation was assessed channel-by-channel for each task (i.e., NTR, CR, P1, P2 and P3) with respect to the baseline condition (i.e., FIX) through a Wilcoxon's test. The FDR method was adopted to correct p-values for multiple testing.
Moreover, statistical difference in the Theta P among the five conditions considered (i.e., NTR, CR, P1, P2 and P3) was also investigated by means of a repeated measures Friedman's test, for each brain area listed in TABLE I. The same test was performed to compare graph measures in five conditions (i.e., FIX, NTR, CR, P1, P2 and P3) for each brain area.
When opportune, Bonferroni's correction was applied to correct the post-hoc analysis for multiple comparisons (k = 6).
Results were considered significant when p<0.05.

A. Spectral Analysis
The frequency analysis was specifically focused on the Theta band power. Fig. 3A depicts the median topographical maps of the relative Theta power during the five analyzed conditions (i.e., NTR, CR, P1, P2 and P3) with respect to the segment of fixation (baseline). These maps highlight a Theta power increase in the frontal and parieto-occipital regions during both the reading (NTR, CR) and the programming (P1-3) tasks with respect to the resting condition (FIX). Moreover, it can be also observed a Theta power increase in the frontal and parietal regions during programming with respect to NTR and CR.
Channel-by-channel statistical analysis revealed a significant increase in Theta power with respect to FIX during P1 on channels belonging to the frontal right, central, temporal left, parietal and occipital regions, and during P2 on channels in the frontal right and occipital regions (Fig. 3B).
Moreover, the Friedman's test applied on the P on a region base, revealed a significant difference among the tasks in the FR (χ 2 = 14, 18, p <0.01) and FM (χ 2 = 12, 44, p <0.05) regions. Specifically, the post-hoc multiple comparison analysis detected a significantly higher increase of Theta P during P1 with respect to CR in the FR region (p=0.03). In the FM region, instead, a theta increased significantly more during P1 with respect to NTR (p=0.04).

B. PSI and Net Flux
The median maps of the net flux measured in the Theta band are depicted in Fig. 3C.
Differently from what happens during the resting phase (i.e., FIX), during NTR, CR, P1, P2 and P3 a clear net-flux pattern can be observed: the parieto-occipital areas assume a source role while frontal areas behave as a sink of the network. Qualitatively, net flux patterns are very similar for all the task-related windows.
Since the net-flux provides an averaged information about how each node behaves in the network, it is not possible to infer that the direction of the information flux goes directly from the parieto-occipital regions to the frontal ones from the observation of the topographical maps.
Such information can be extracted by means of the median behavior of the PSI index depicted in Fig. 4. This representation shows that parietal nodes send information both towards frontal and occipital nodes during all tasks. Moreover, occipital nodes, in addition to receiving information from the parietal area, also send information towards the frontal nodes.
C. Network and Graph Analysis 1) Reduced Networks: The flux direction was also analyzed by means of the RNs (Paragraph II.F.I) showed in Fig. 5. The hypothesis that the information flows from the parietal to the frontal areas was confirmed by their median trends.
The presence of a significant difference in the median connectivity between the parietal region and the frontal ones (i.e., P→F) among the five conditions was revealed by a Friedman's test (χ 2 = 19, 86, p<0.05). A successive post-hoc test for multiple comparison with Bonferroni's correction revealed a significant increase of P->F between the FIX segment and P1 (p=0.01).
2) Betweenness of Centrality: The BC of the FM region decreased passing from the resting condition (FIX) to the cognitive tasks (i.e., NTR, CR, P1, P2 and P3). This decrease was particularly enhanced during the programming windows. In the P region, instead, an increased BC was observed during programming with respect to the other tasks (Fig. 6).
Friedman's test revealed significant differences of BC among the protocol conditions only in region FM (χ 2 = 33.5, p≪0.001). Concerning the FM area, post-hoc comparisons detected a significant decrease of BC from the FIX segment to P1, P2 and P3 (p = 0.02, p = 9.9e-4 and p = 9.9e-4, respectively).
3) In and Out Degree: ID trends showed a general increase in the frontal and central areas during reading and   programming with respect to FIX, while OD mainly decreased during programming in the parietal region. Specifically, statistical analysis revealed significant differences for the ID in the regions FM (χ 2 = 16.8, p<0.01), and TR (χ 2 = 18.39, p<0.01). A significant increase was detected in FM during P1 with respect to FIX (p = 0.003), despite NTR and CR showed a trend comparable to P1 (Fig. 7A). Moreover, ID increased significantly in P1 and P3 with respect to FIX (p = 0.01 and p = 0.02, respectively. Significant results in the FR (χ 2 = 19.64, p<0.01) and P (χ 2 = 15.06, p<0.01) areas were obtained for the OD.
In region FR, OD increases significantly during P1, P2 and P3 with respect to FIX (p = 0.02, p = 0.01, p = 0.009). In the P area a significant (p = 0.01) increase in P1 with respect to FIX was detected (Fig. 7B). 4) Clustering Coefficient: No significant differences were observed in terms of CC among the five conditions considered. 5) Global Efficiency: Neither significant nor qualitative differences were observed in terms of GE among the five conditions considered.

IV. DISCUSSION
The main goal of the present study was to shed lights on the brain mechanisms that support software programming by investigating the modulation in power and functional connectivity patterns of cortical Theta rhythm. Indeed, Theta oscillations have been shown to support complex cognitive processing and, in accordance with previous studies Fig. 7. In Degree distributions for the FM brain region (A) and Out Degree distributions for the P brain region (B) computed during FIX, NTR, CR, P1, P2 and P3. Asterisks indicate significant difference between tasks.
[17], [18], the contribution of the Theta rhythm to the total EEG power increased both during reading tasks (NTR and CR) and during programming tasks (P1, P2 and P3), while for the first time, we also described changes in network topology, suggesting a task-related adaptation of the supporting functional system.

A. Modulation of Theta Spectral Power and Network Arrangement
Theta power increase was particularly enhanced in the frontal and posterior areas during coding. These results suggest the involvement of WM processes during programming tasks.
Theta activity, indeed, seems to support WM, which is a brain system involved in the manipulation of verbal and visual inputs temporary stored in the short-term memory [34]. WM processes are involved in several cognitive tasks, such as language processing and mathematical calculation. Therefore, the dominant role of Theta oscillations, observed in this study, may confirm the involvement of such processes also in software programming. Moreover, the major enhancement of Theta power that occurred during programming, with respect to the reading segments, could be associated to an increased task demand and allocation of WM resources [35]. Software programming, indeed, is a very complex task, which foresees both language comprehension, logic and symbol manipulation processes. In line with the spectral analysis results, the PSI, computed in the Theta band, highlighted the main role of frontal and posterior regions in cognitive tasks. WM processes, indeed, result from the dynamic interplay of specific and distributed fronto-parietal brain regions, responsible for information encoding, data retrieval and sensory input processing [34]. As suggested by Ratcliffe and colleagues, frontal regions are supposed to be involved in executive functions, while the posterior ones seem to be responsible for information processing and maintenance [35]. Moreover, a modulation of parietal-frontal integration in the Theta band was observed during a letter-based sustained attention task [36]. Indeed, our results showed that the preferential directionality of the information flux goes from posterior brain regions to frontal ones (Fig. 5).
Specifically, parietal nodes seem to assume a prominent role in the Theta network, sending information both to frontal and to occipital regions. Interestingly, as shown by the RNs, during the first and final windows of the coding assignment (i.e., P1 and P3) the strength of the normalized connection from parietal electrodes to frontal ones increased with respect to all the other considered segments of the protocol (i.e., FIX, NTR, CR and P2). This behavior could be again associated to an increase in cognitive demand. Indeed, programming involves tasks like NTR and CR, but also problem-solving and code generation. For this reason, the amount of information to be processed and the need of memory retrieval may be magnified. Interestingly, the increase of PSI links directed from the parietal region to frontal ones was significant only between FIX and P1. In fact, during the first window of coding all subjects were reading instructions and retrieving their knowledge to comprehend and solve the assignment. Therefore, the higher degree of parietal to frontal connection observed during P1 could be due to the occurrence of both language comprehension and exercise planning. Though spectral analysis and PSI results showed that both reading segments and programming involve Theta brain fronto-parietal networks, graph analysis enhanced some task and time dependent differences in their topology.
The most interesting result is associated to the measure of centrality, which decreased in frontal regions and increased at the parietal ones during coding. This result is explainable with an increase in the number of network paths that go through the parietal nodes and a contemporary drop of the ones passing through the frontal areas.
Moreover, while the number of frontal inward connections increased for all tasks with respect to FIX, the number of outward connections from the parietal region showed a pronounced increase during the coding windows. Finally, a significant negative correlation between the P→F and the perceived mental effort was found (r = −0.5, p = 0.01).
This result confirms the importance of the parietal region in programming tasks. More in detail, the higher is the information flow that starts from parietal nodes and arrives to the frontal one, the lower is the effort perceived by the subject, suggesting an association between parietal source activity and the performance of the network. Our hypothesis is that pre-existing networks historically used for older cognitive processes specialized over years to accomplish programming tasks, such as programming and code comprehension. We suggest that this specialization manifests in brain connectivity dynamics, that is the ability of brain networks to reorganize their topology in function of the executed task.

B. Limitations of the Study
The small sample size considered (i.e., 11 subjects) obviously is the main limitation of the present research, despite this, we obtained promising results. To obtain more robust results and to further investigate the described findings, the number of involved programmers should be increased.
The protocol we adopted, if, in one hand, has the merit of reproducing conditions similar to a real working scenario (event-free programming task), the use of sensors for biosignals acquisition, in the other hand, reduced the comfort of the participants, potentially introducing physical annoyance. despite their noninvasiveness.
Additionally, this study did not involve a control group (i.e., subjects not expert in programming), but only a group of experienced programmers who were tested during different conditions, while the baseline was the "control" state. For this reason, it was not possible to perform the assessment of proficiency-dependent modulations of brain connectivity. This is a very interesting aspect which requires the comparison of groups homogeneous for their proficiency and could be object of a future study.
Finally, as concern the methodological approach we adopted in the signal processing, we should point out the particular way the Phase Slope index computes the directional information. Indeed, the PSI is an antisymmetric index, thus it identifies only one direction of the analyzed connection between two nodes. In this way, PSI does not describe the existence of bidirectional couplings, which, conversely, indexes based on the GC can identify (GC xy ̸ = GC yx ). Nevertheless, PSI is the index that provides a total and net information about directionality and was here selected for its straightforward interpretability and for its robustness against VC effect.

V. CONCLUSION
The current study presents a step forward in the understanding of the neural mechanisms supporting software programming tasks in terms of brain rhythms and interactions. Specifically, our results revealed that tasks of natural text reading, code reading, and realistic programming all recruit Theta fronto-parietal networks. This evidence could be reasonably associated to the involvement of working memory processes. Graph analysis showed that some task dependent changes occur in the topology of involved networks. The main finding resides in the ability of the networks to dynamically adapt to task changes and requirements. Specifically, this network reorganization mainly involved parietal areas during programming, which increased their connectivity in terms of centrality and degree in the network. Moreover, a negative correlation between the parietal out degree and the mental effort perceived by the subjects was observed. These results suggest that the topology of the network in the Theta band, especially at the parietal level, is promising in the discrimination of the performed task and of the subjects' cognitive effort. Since Theta brain networks seem to modulate their topology in function of the task requirements, further knowledge in this field could be reached in the future by a continuous monitoring of the dynamic interplay of fronto-parietal regions, with the aim of studying event-related changes in network topology. This dynamicity of brain networks can be reliably tracked by means of EEG signals, given their high temporal resolution [37].
Moreover, further knowledge in the field could be reached by conducting the same study on groups of programmers with different levels of proficiency, in order to highlight possible differences in Theta functional networks organization due to experience.