{'fraction_non_alphanumeric': 0.03038454089647307, 'fraction_numerical': 0.026360525919390587, 'mean_word_length': 4.58356361888742, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Using Machine Learning to Predict the Evolution of Physics Research

The advancement of science as outlined by Popper and Kuhn is largely qualitative, but with bibliometric data it is possible and desirable to develop a quantitative picture of scientific progress. Furthermore it is also important to allocate finite resources to research topics that have growth potential, to accelerate the process from scientific breakthroughs to technological innovations. In this paper, we address this problem of quantitative knowledge evolution by analysing the APS publication data set from 1981 to 2010. We build the bibliographic coupling and co-citation networks, use the Louvain method to detect topical clusters (TCs) in each year, measure the similarity of TCs in consecutive years, and visualize the results as alluvial diagrams. Having the predictive features describing a given TC and its known evolution in the next year, we can train a machine learning model to predict future changes of TCs, i.e., their continuing, dissolving, merging and splitting. We found the number of papers from certain journals, the degree, closeness, and betweenness to be the most predictive features. Additionally, betweenness increases significantly for merging events, and decreases significantly for splitting events. Our results represent a first step from a descriptive understanding of the Science of Science (SciSci), towards one that is ultimately prescriptive.


Introduction
We all become scientists because we want to create an impact and make a difference, to the lives of those around us, and also to the many generations that are to come. We all strive to make choices in the problems we study, but not all choices lead to breakthroughs. There is actually a lot more about scientific breakthroughs that we can try to understand. For one, science is an ecosystem of scholars, ideas, and papers published. In this ecosystem, scientists can form strongly-interacting groups over a particular period to solve specific problems, but later drift apart as their interests diverge, or due to the availability or paucity of funds, or other factors. The evolution of these problem-driven groups is more or less completely documented by the papers published as outcomes of their research. By analysing groups of closely related papers, researchers could extract rich information about knowledge processes [1][2][3] . The potential to map scientific progress using publication data has attracted enormous interest recently [4][5][6] . However, compared to the study on science at the level of individual papers [7][8][9] and at the level of the whole citation network [10][11][12] , where a lot of work has already been done, the research on science at the community level is still limited 1, 3 .
In a recent paper, Liu et al. demonstrated the utility of visualizing and analysing scientific knowledge evolution for physics at the aggregated mesoscale through the use of alluvial diagrams 3 . In this picture, papers are clustered into groups (or communities) and these groups can grow or shrink, merge or split, new groups may arise while the others may dissolve. This shares a very strong parallel with what some researchers discovered in social group dynamics 13 . More importantly, many breakthroughs were made by scientists absorbing knowledge from other fields, often in a very short time. On the alluvial diagrams, these knowledge transformations manifest themselves as merging and splitting events. Clearly, funding agencies, universities and research institutes would want to promote growing research fields, and particularly those where breakthroughs are imminent. This is why it is important to be able to predict the future events. Liu et al. 3 attempted this in their paper, by analysing the correlation between event types and several network metrics. Unfortunately, such predictions are very noisy.
(BCN) reflects the relation between present publications while the co-citation network (CN) represent the relation between papers which have strong influence on recent publications. In this way, we can detect communities over the years, and study how they evolve year by year, see the Methods section for details on BCN and CN.
After building BCN and CN, the Louvain method was used to extract the community structures. By checking the Physics and Astronomy Classification Scheme (PACS) numbers of the papers in these communities, we have shown that the BCN communities are meaningful and reflect the real structure of the scientific communities 3 . We found indeed that papers in the same community are really focus on the closely related topics. For the CN communities, this validation is tricky because of two problems: (i) the old Physics Review papers have no PACS numbers, and (ii) PACS was revised several times, so the same numbers in different versions can potentially refer to different topics, or the same topics are referred to by different numbers in different versions. Nevertheless, systematic validation seems to be impossible although a quick check on some CN communities after 2010 suggests that CN community structure also reliably reflects the actual scientific community. We refer to these validated units of knowledge evolution as topical clusters (TCs) in this paper.
In Figure 1, we provide the alluvial diagram that depicts the evolution of TCs within the BCNs for the period between 1981 to 2010. The equivalent alluvial diagram for the CNs is shown in Figure S2 in Supplementary Information (SI). In both alluvial diagrams, we visualized the sequences of TCs, their inheritance relations, which can be intimacy indices (for the BCN communities), fraction of common members or inclusion measures (for the CN communities), and the evolution processes they undergo, see the Methods section for more details. The events (changes) that we can discern from the alluvial diagram (shown in Figure 1) are analogous to those recognized in social group evolution 13 . They represent forming, dissolving, growing, shrinking, merging and splitting. We found in Liu et al. that the prediction of such events is hard, since the correlation between them is nonlinear and complex. This challenge is addressed in the following section by tapping into the power of machine learning.

Event labelling
The GED method takes into account the size and the similarity between groups (TCs) in the consecutive time frames in order to label groups change (assign event type). There are four events considered in this work: • continuing-a research field is said to be continuing when the problems identified and solutions obtained from one year to another are of an incremental nature. It is likely correspond to the repeated hypothesis testing picture of the progress of science proposed by Karl Popper 24 . Therefore, in the CN, this would appear as a group of papers that are repeatedly together cited year by year. In the BCN, this shows up as groups of articles from successive years sharing more or less the same reference list.
• dissolving-a research field is thought to disappear in the following year if the problems are solved or abandoned, and no new significant work is done after this. For the CN, we will find a group of papers that are cited up to a given year, but receives very few new citations afterwards. In the BCN, no new relevant papers are published in the field, hence, the reference chain terminates.
• splitting--a research field splits in the following year, when the community of scientists who used to work on the same problems, start to form two or more sub-communities, which are more and more distant from one another. In terms of the CN, we will find a group of papers that are almost always cited together up till a given year, breaking up into smaller and disjoint groups of papers that are cited together in the next year. In the BCN, we will find the transition between new papers citing a group of older papers to new papers citing only a part of this reference group.
• merging--multiple research fields are considered to have merged in the following year when the previously disjoint communities of scientists found mutual interest in each other's field so that they solve the problems in their own domain using methods from another domain. In the CN, we find previously distinct groups of papers that are cited together by papers published after a given year. In the BCN, newly published papers will form a group commonly citing several previously disjoint groups of older papers.
The GED method has two main parameters (alpha and beta), which are the levels of inclusion that groups in the consecutive years have to cross in order to be considered as matching groups. We have applied the GED method with the wide range of these parameters from 5% to 100%. The characteristics of the considered networks required us to set the alpha and beta thresholds to very low values, i.e. 30% for the BCN, and 10% for the CN, see SI for more details. In total, we have obtained 479 various events for the BCN, and 492 events for the CN, which are our observations and the labels in the prediction part of our study. In both networks, the events distribution was imbalanced with the continuing event dominating over all other types, see Figure 2A1 and Figure 2B1. Each block in a column represents a TC and the height of the block is proportional to the number of papers in the TC. For clarity reason only TCs comprising more than 100 papers are shown. TCs in successive years are connected by streams whose widths at the left and right ends are proportional to the forward and backward intimacy indices. The colours inside a TC represent the relative contributions from different journals.

Future events prediction
The machine learning approach to prediction requires dividing the data into two parts: the training data set and test data set. The training data is used to learn classifier, which can then label events in the test data. The labelled values are compared with the event labels and the prediction performance is calculated. More than 450 observations were used to train the classifiers. Each observation contained 77 features (preselected from the initial 100) divided into three categories: microscopic features (related to nodes in the group, e.g. node degree), mesoscopic features (related to the entire group, e.g. the group size), and macroscopic features (related to the whole network, e.g. network density). Mesoscopic features calculated for individual nodes are commonly aggregated for all nodes from the group, e.g. average node degree or betweenness in the group. See the SI for the complete list of features used. To automatically select the best classification algorithm (model) as well as its hyper-parameter settings to maximize the prediction performance, the Auto-WEKA software package 25 was utilized. For each network, we ran the Auto-WEKA for 48 hours, which allowed us to validate nearly 20,000 configurations per network. The metric being maximized was the F-measure, commonly used for multi-class classification. The overall classification quality was calculated as the average F-measure for all event types, treating them as equally important.
The predicted output variable (event labels) had an imbalanced distribution. Commonly, classifiers tend to focus on the dominant event type (class), which is very well predicted, but at the expense of the minority event types. For the imbalanced BCN data set, the best performance was achieved with the Attribute Selected Classifier (with the SMO as base classifier), which performs feature selection 26 . The percentage of the correctly classified instances was 80.6%, while the average F-measure was only 0.50 due to classifier focusing on continuing, which was the most frequently occurring event type, see Figure 2A. For this event, the F-measure value was equal to 0.89, and only 7 events out of 352 were incorrectly classified. The worst classified was the splitting event, whose F-measure was only as of 0.11. Most of the splitting events were incorrectly classified as continuing (31 out of 33 events). The second worst was merging, with F-measure 0.35. Again, the majority of the merging events were wrongly classified as continuing events: 38 out of 56. Interestingly, the splitting and merging events were never cross-classified mistakenly. For the imbalanced CN data set, the best performance was achieved with a lazy classifier, which uses locally weighted learning 27 . The percentage of the correctly classified instances was 73.3%, while the average F-measure was only 0.53, again due to the classifier concentrating on the dominating continuing event type, see Figure 2B. The F-measure value for the continuing event was only 0.83, however, as many as 50 continuing events (out of 337) were wrongly classified as dissolving. Alike to BCN, many splitting and merging events were incorrectly classified as continuing: 17 out of 22 events, and 24 out of 46 events, with F-measure equal to 0.30 and 0.42, respectively.
By balancing the imbalanced training data sets (i.e. by equally sampling them), we force the classifiers to pay more attention to the features rather than to the number of occurrences of the particular majority event type. As a result of balancing data sets, the previously minor event types (dissolving, merging, and splitting) were predicted much better, but with a significant drop in performance of the continuing event classification. More importantly, by balancing the data sets we increased the overall prediction quality by over 20%. For the balanced BCN data set, the best performance was achieved by means of the boosting-based classifier AdaBoost with the Bayes Net as the base model. The percentage of the correctly classified instances was 62.0% and the average F-measure was 0.61. The biggest sources of errors were merging events, which were wrongly classified as continuing and dissolving, as well as continuing wrongly classified as splitting. The best classified event was dissolving (only 4 mistakes in 27 classifications, the overall score 0.79) followed by the splitting event (6 mistakes in 27 classifications, overall F-measure 0.70). For the balanced CN data set, the Attribute Selected Classifier (with the PART as base classifier) provided the best results-the percentage of the correctly classified instances was 69.32%, while the average F-measure was 0.69 28 . The dissolving, merging, and splitting events were classified very well with the F-measure values equal to 0.79, 0.82, and 0.75 respectively. Most of the continuing events were wrongly classified as splitting (13 out of 22), which resulted in lower F-measure value 0.40.
What is interesting for us to note is that the prediction results for the CN being slightly better than for the BCN. A possible explanation is that for the CN we used a richer similarity measure containing users importance information. Thus the event tracking and therefore the ground truth could be more accurate. Overall, the prediction quality expressed by the average F-measure was very good for the imbalanced as well as for the balanced data sets, as the baseline results obtained with the ZeroR classifier were much worse: F-measure 0.21 for both, BCN and CN, imbalanced data sets, 0.18 for the balanced BCN and 0.14 for the balanced CN. For each data set different classifier turned out to be the best, however most models were wrapped with the boosting or meta classifiers.

Predictive feature ranking
The feature selection technique is used in machine learning to find the most informative features, to avoid classifier overfitting, to eliminate (or at least to reduce) the noise in the data as well as to provide some explanations about phenomena 29 . By repeating the feature selection 1000 times, we obtained 1000 sets of selected features. Next, we calculated how many times each feature has been selected, thus, receiving the ranking of the most often selected features. For the BCN, the context-based The context-based features (number of papers published in a given journal) turned out to be the most informative, followed by the microscopic structural measures, especially closeness, degree and betweenness.
features dominated the ranking. It referred, especially the number of papers from the Physical Review E, Physical Review B and Physical Review A, see Figure 3A. Beside the context, the network features based on degree, betweenness, size and closeness measures were most informative, which tells us that the structural properties are as important as context awareness. The context-based feature, i.e. the number of papers published in Review of Modern Physics, was the most often selected for the CN data set. It is followed by closeness-and degree-based features in the ranking, see Figure 3B. For both networks macroscopic features were ranked rather low, which suggests that the general network profile is not very important, perhaps because of the smooth changes in the entire network. Surprisingly, the dynamic features, e.g. related to the average age of references (for BCN) and age of articles (for CN) did not show informative value and were ranked very low for both networks. The rankings were validated in the additional two years of data available (2010-2012). The prediction was performed twice: (i) using all features, and (ii) using the top 10 ranked features only. Selecting only the top 10 features, boosted the quality of the prediction by 11% for the CN, and by 2% for the BCN, which underlines the necessity of the feature selection process.

Changes to the Betweenness Distributions Associated with Merging and Splitting Events in BCN
Having the list of best predictive features, Figure 3, we can analyse some of them more in-depth to look for early warning signals. Basically, we believe that scientific knowledge evolves slowly, and this slow evolution drives the evolution of citation patterns. Therefore, there must be specific changes in citation patterns that precede merging and splitting events. Besides the number of PRE papers in a TC, the sum network betweenness is also a strongly predictive feature, see Figure 3A. This suggests that we should look at the betweenness of papers in the BCN more carefully. The betweenness of the node denotes what percentage of shortest paths between all pairs of nodes in the network passes a given node. Values of nodes' betweenness can be aggregated (sum, average, max, min) for all nodes in the TC, as what we list in Table S1. However, in this section we only focus on the distribution of original node betweenness. Naively, when we consider the part of the BCN adjacency matrix corresponding to two TCs that ultimately merged, we expect to find few links between TCs at first. But as the number of links between TCs increase over time, the modularity-maximizing Louvain method will eventually merge the two TCs into a single TC. This is shown schematically in Figure 4, where in general betweenness will increase on average with time as the two TCs merge. In reality, there are always links between TCs, and the numbers and strengths of these links fluctuate over time. To develop a more quantitative description of the merging events outlined in Figure 1, as well as splitting and continuing events, we focus on five events going from 1999 to 2000, shown in Table 1.  In Table 2, we show the 25th, 50th and 75th percentiles of the papers in these smaller groups, compared to those of 1849 papers in 1999.01 and 344 papers in 1999.02. As we can see, the 25th, 50th, 75th percentile betweenness in the connecting parts (1999.01a and 1999.02a) are all higher than the 25th, 50th, 75th percentile betweenness in the non-connecting parts (1999.01b and 1999.02b). More importantly, these percentile betweenness are higher than the 25th, 50th, 75th percentile betweenness of the TCs 1999.01 and 1999.02 themselves. To test how significant these quartiles are in 1999.01a, we randomly sampled 164 betweenness values from 1999.01 10 6 times, and measured the quartiles of these samples. When we draw random samples from a TC, the 25th percentile, the 50th percentile, and the 75th percentile, depends on the size of the TC. There is more variability in these quartiles in smaller samples than they are in larger samples. Therefore, in the test for statistical significance, the observed quartile has to be tested against different null model quartiles for samples of different sizes. To do this, we draw samples with a range of sizes from the same set of betweenness, and for a given quartile (25%, 50%, or 75%), fit the minimum quartile value against sample size to a cubic spline, and the maximum quartile value against sample size to a different cubic spline. With these two cubic splines, we can then check whether the observed quartile value for a sample of size n is more than or less than the null model minimum or maximum using cubic spline interpolation. From the histograms shown in Figure 6(A), we see that the betweenness quartiles of 1999.01a are statistically larger than random samples of the same size from 1999.01, at the level of p < 10 −6 , which means the papers in 1999.01a have significantly larger betweenness than other papers in 1999.01.
We also checked the statistical significance of the larger betweenness values of 1999.02a, against 10 6 random samples of the same length (144) from 1999.02. From Figure 6(B), we can derive that the quartiles of 1999.02a are only a little larger than the tails of the quartile histograms of the random samples, but their statistical significance is still at the level of p < 10 −6 . When a TC splits into two in the next year, we expect the links between two parts a and b in the TC to have thinned out to the point that the modularity Q of the whole is lower than the modularities Q a and Q b of the two parts. However, in general, we would not know how to separate the TC into the two parts a and b. Fortunately, for the 1999.01 → 2000.02 + 2000.03 splitting event, we also know the part 1999.01a, which merged with 1999.02a, became 2000.03. Therefore, we might naively expect 1999.01b to be the part that split from 1999.01 to become 2000.02. If we test the quartiles of 1999.01b, against random samples of the same size from 1999.01, we find the histograms shown in Figure 6(C). As we can see, the betweenness quartiles of 1999.01b are quite a bit lower than the typical values in 1999.01, but this difference is statistically not as significant as the quartiles of 1999.01a. Thinking about this problem more deeply, we realized that while papers in 1999.01b have no references in common with 1999.02, some of them do share common references with 1999.01a. Let us call these sets of papers 1999.01aα (papers do not have references in common with papers in 1999.01b), 1999.01aβ (papers have references in common with papers in 1999.01b), 1999.01bα(papers have references in common with papers in 1999.01a), and 1999.01bβ (papers that do not have references in common with papers in 1999.01a). In Figure 6(D), we learn from the histograms that the betweenness quartiles of 1999.01bα are indistinguishable with random samples of the same size from 1999.01. On the other hand, from the histograms in Figure 6(E), we find out that while the lower betweenness quartile of 1999.01bβ is indistinguishable with the random samples of the same size from 1999.01, its median and upper quartile are both on the low sides of the random sample distributions. This suggests a split of 1999.01 to (1999.01a + 1999.01bα) + 1999.01bβ .
Just to be safe, we also checked the betweenness quartiles of 1999.01aα and 1999.01aβ , against random samples of the same sizes from 1999.01. As we can see from Figure 6(F) and (G), the lower quartiles and medians are lower than those obtained from random samples, but the upper quartiles are decidedly higher. However, the difference between 1999.01aα and 1999.01aβ is not as obvious as difference between 1999.01bα and 1999.01bβ , one possible reason is the smaller sample size (17, 147 vs. 907, 778). Again, these results are consistent with the picture that the rise in betweenness in 1999.01a is driving  Although a small part split off from each of 1999.11 and 1999.12, the main event associated with the two TCs was a symmetric merging. Looking again into the relevant parts of the BCN, we found 299 out of 1014 papers in 1999.11 coupled to 347 out of 988 papers in 1999.12, and we call them 1999.11a and 1999.12a, respectively. As we can see from the histograms in Figure 6(H) and (J), the betweenness quartiles in 1999.11a and 1999.12a are significantly higher than one would expect from random samples of 1999.11 and 1999.12. Simultaneously, the betweenness quartiles in 1999.11b and 1999.12b are significantly lower than in random samples of 1999.11 and 1999.12 (see Figure 6(I) and (K)). Therefore, what we are seeing here might be the early warning signals of merging, as well as that of the asymmetric splitting. So far we have learnt that a decrease in betweenness within a TC signals a possible split, whereas an increase in betweenness of the part of the TC coupled to another TC signals a merger between the two TCs. For this story to be consistent, we must not see these signals in the continuing events 1999.04 → 2000.06 and 1999.13 → 2000.16. However, if we go through the full BCN, we find that 370 out of 389 papers in 1999.04 and 308 out of 319 papers in 1999.13 are coupled to papers outside of these TCs, which suggests the possibility of merging or splitting. However, as we can conclude from Table 3, while the lower betweenness quartiles of the coupling parts of 1999.04 and 1999.13 with other TCs may be significantly larger than those of random samples of the two TCs, the highest betweenness quartiles are never significantly larger. Therefore, at the same level of confidence that we have set for the precursors of merging between 1999.01 and 1999.02, as well as between 1999.11 and 1999.12, we have to say that there is no significant precursors for 1999.04 and 1999.13 to merge with other TCs. What about splitting then? A TC is likely to split into two if at least one of two parts has reduced betweenness. We see in Table 3 that betweenness in the coupling parts of 1999.04 and 1999.13 are not significantly lower than those of random samples. Therefore, we look at the non-coupling part, i.e. papers in 1999.04 and 1999.13 which have no references in common with papers in other TCs, but they may have common references with papers in the same TCs. We call these non-coupling parts 1999.04b and 1999.13b, respectively (the bottom row in Table 3). Only the top betweenness quartile of 1999.04b falls below that of random samples from 1999.04 in Table 3. Therefore, the early warning for a splitting event in the next year is not strong enough. For 1999.13b, on the other hand, all three betweenness quartiles fall below that of random samples from 1999.13, even after we have accounted for the small size of 1999.13b. This suggests that the probabiiity of a splitting event next year is high, but 1999.13 continued on to 2000.16, which thereafter continued to 2001 without merging or splitting. This might be because additional conditions, like the size of TC being large, must be satisfied before a splitting can occur.

Discussion
During the past two decades, researchers have made a lot of efforts to understand the system of science. Many problems are solved, however the understanding of interactions between different fields is still limited. Investigating the temporal network (BCN, CN) and their community structures, we are able to measure and quantify the complex interaction between different fields particularly in physics over time. Naturally, we would like to have a predictive power based on this picture. However, the correlation between network structure and evolution events is nonlinear and complex. Therefore we turn to machine learning techniques, which have shown a great power to solve predictive problems that are hardly to be solved using traditional statistical methods. To our knowledge, this is the first study that utilizes both machine learning and network science approaches to predict the future of science at the community level.
To be able to identify changes in TCs we needed to define time windows used for network creation and community detection. The natural choice for bibliographical data was the usage of single years, since the publishing process may last many months. Obviously, another granularity may be considered like multiple years, e.g. 2 or 5 years. In our approach, i.e. both for BCN and CN, every citation has the same importance. However, there are some other concepts like fractional counting of citations 30 . It assumes that the impact of each citation is proportionate to the number of references in the citing document. Additionally, it can be differentiated depending on e.g. the quality of the journal. For the CN we have calculated the similarity between groups in the consecutive time windows in two ways: (i) using the plain relative overlap measure, and (ii) using the inclusion measure based on Social Position. The idea was to enrich evolution data with the structural information occurring between the nodes. It turned out that both measures provided the similar labelling, but the evolution tracking with the Social Position information produced slightly better initial prediction. Therefore, the study was continued only for the inclusion measure,see SI for more information.
We decided to analyse more in-depth only on one feature describing structural profile of TCs, namely node betweenness. It was primary caused by the limited amount of resources and complexity of analyses. The entire process required much human assistance and could not have been easily automated. In our experiments, we utilized the raw, imbalanced or artificially flattened-balanced data sets. However, depending on the study purpose, we can bias some classes we are more interested in e.g. split. It can be achieved either by means of appropriate balancing-sampling for the learning set, or reformulating the problem into the binary question-is split expected (true) or not (false). As of now, the betweenness analysis is still limited to several case studies, in future a more rigorous framework will be desired. The idea of analysing science by discovery of knowledge changes is general and can be applied to all bibliographical data containing citations. We focus solely on APS journals, however, also papers indexed by PubMed, Web of Science or Google Scholar may be studied.

Methods
The entire analytical process consists of several steps that are primary defined by the Group Evolution Prediction (GEP) framework. First, the bibliographic coupling network (BCN) and co-citation network (CN) are extracted from the references placed in the papers from a given time window, see Figure 7, and this is carried out separately for each period. As a result, we get a time series of BCNs/CNs. Next, paper groups called topical clusters (TCs) are extracted using the Louvain clustering methods, independently for each BCN/CN in the time series. Each group is described by the set of predictive features. Having TCs for consecutive periods, we were able to identify changes in TC evolution using the Group Evolution Discovery (GED) method that appropriately labels the TC changes, see below.
Independently, the features ranking and its validation were performed to find the most valuable TC measures. Based on this ranking, a structural measure node betweenness was selected for the more in-depth studies as the early signal for splitting or merging.

GEP method
The Group Evolution Prediction (GEP) method is the first generic approach for the prediction of the evolution of groups 16 , in our case groups correspond to TCs. The GEP process consists of six main steps: (1) time window definition, (2) temporal network creation, (3) group detection, (4) group evolution tracking, (5) evolution chain identification and feature calculation, and (6) classification using machine learning techniques. Thanks to its adaptable character, we were able to apply it to the BCN and CN differently. For the group (TC) detection in both networks, we applied the Louvain method 31 . The group evolution tracking was performed with the GED method (see below), but we used different similarity measures for each network BCN and CN, see below. The set of features describing the group at a given time windows was adjusted to our networks, as some of the features defined in the GEP method were not applicable in our case. We also introduced some new, dedicated measures appropriate for bibliographical data, see SI for the complete list. Finally, we applied the Auto-WEKA tool to find the best predictive model and its parameters from the wide range of all possible solutions. The commonly known average F-measure was used as a prediction performance measure. Figure 7. The process of building a Bibliographical Coupling Network (BCN) and Co-citation Network (CN) from the citation bipartite network for a given period-year t. Both BCN and CN are undirected and weighted; the weights denote the number of shared citations (BCN) or co-citing papers (CN). Separate topical clusters are extracted for BCN (C 1 , C 2 ) and CN (C 3 , C 4 ). Nodes with numbers are papers from a given period being considered and nodes with letters are their references.

Bibliographic coupling network (BCN) and co-citation network (CN)
In the BCN and CN, nodes represent papers and undirected but weighted edges denote the bibliographic coupling strengths and co-citation strengths, respectively. That is, if two papers share w common references, the BCN edge between them would have a weight of w. For example, papers 1 and 2 in Figure 7 share three citations: A, B, and C, whereas papers 3 and 4 commonly cite only one paper-E. On the other hand, if two papers are cited together by w papers, the edge between them in the CN receives weight w . Papers A and B are cited together by two other papers: 1 and 2, but papers B and C by three, i.e. additionally by paper 3. Both BCN and CN are temporal networks, in which the nodes are all papers published within a specific time window (BCN) or papers cited within a given time window (CN). We assume that the reasonable time window for bibliographical data is one year to facilitate the analysis of changes in scientific knowledge, i.e. changes in topical clusters year by year. For the BCN, only the giant component, which in most cases occupies 99% of the whole BCN, will be considered for the TC detection and evolution analysis. For the CN, we do not use all papers cited in the given time window because most of them are cited only a small number of times, and thus they have little influence on the broader knowledge evolution. Therefore, we rank all available N papers p 1 , p 2 , . . . , p N in the descending order by the number of times they are cited in this time window (year): Next, we choose the top n papers p 1 , p 2 , . . . , p n , that totally gathered 1 4 of all citations, i.e. such that n < N is the smallest integer to satisfy ∑ n i=1 f i ≥ 1 4 ∑ N j=1 f j . The data we used in this paper is the APS data set, consisting of about half a million publications between 1893 to 2013 and six million citation relations among them 32 .

Intimacy indices
To analyse the evolution of TCs, we need to match them from the consecutive years. The set of cited papers to large extent overlaps year by year, so for the CN, we can use the regular approach proposed together with the GED method, see below and Brodka et al. 20 . For BCN, however, there is no overlap at all between papers published in the successive years because every paper can be published only once and in only one year. Even if we do not have the corresponding papers in TCs from two BCNs, i.e. two years, the papers' references overlap each another. Therefore, we can measure the similarity of their reference pools to reflect their inheritance. For that purpose, we introduced the forward intimacy index and backward intimacy index in Liu et al. 3 . The idea behind intimacy indices is that the references related to a particular topic change gradually. The forward intimacy

14/24
index I f mn and the backward intimacy index I b mn between TCs C t m in year t and C t+1 n in year t + 1 are defined in following: , (1) Here the TCs at t and t +1 are C t = {C t 1 , ...,C t m , ...,C t u } and C t+1 = C t+1 1 , ...,C t+1 n , ...,C t+1 v , and we denote the references cited by papers in

GED method
The Group Evolution Discovery (GED) method 20 was used for tracking group evolution for historical cases--to learn the classifier and for testing cases to validate classification results. The GED method makes use of the similarity between groups in the following years as well as their sizes to label one of six event types: continuing, dissolving, merging, splitting, growing, shrinking. However, we have adapted the GED method to label only four types of events: continuing, dissolving, merging, splitting, as these are the most important to us. The other two (growing and shrinking) are covered by continuing. In general, the GED method allows us to use various metrics as a similarity measure between groups. Therefore, the intimacy indices defined in Equation 1 were used for the BCN to match similar groups in the consecutive time windows. However, the original GED inclusion measures were used for the CN. It means that the similarity between two groups from two successive time windows is reflected by the inclusion measure, which is calculated for two scenarios: inclusion I(C t n ,C t+1 m ) of a group C t n from time window t in another group C t+1 m from time window t + 1 (forward, Equation 2), and inclusion I(C t+1 m ,C t n ) of this second group C t+1 m from t + 1 in the first group C t n from t (backward, Equation 3). The inclusion measure makes use of the Social Position SP(p), which is a kind of weighted PageRank. It denotes an importance of paper p being cited among all other papers 23 . The inclusions for CN are defined as follows: If both inclusions (CN) or both intimacy indices (BCN) are greater than the percentage thresholds alpha and beta (the only parameters in this method), the method labels the event continuing. If at least one inclusion or one intimacy index exceeds one of the thresholds, the splitting and merging events is considered, the proper event is assigned depending on the number of similar groups in t and t + 1. If both inclusions or both intimacy indexes are below the thresholds, i.e. the group has no corresponding group in the next time window, the dissolving event is assigned.

Feature Ranking
Rankings of the most prominent features was obtained by repeating the feature selection 1000 times using a basic evolutionary algorithm 29 , as proposed in Saganowski et al. 16 . The rankings were received for the 30 years span (1981-2010). Next, only top 10 features were selected to described TCs in two additional years (2010-2012) and predict TC evolution. The results revealed the superiority of feature selection compared to the raw approach with all features engagement.

The alpha and beta thresholds for event labelling
Two groups in the consecutive time windows are considered similar if at least one of their inclusion measures is greater than alpha or beta parameters. In other words, the alpha and beta parameters are thresholds which have to be satisfied to assign an event between two groups. The theoretical range of values for alpha and beta is between 0% and 100%. However, the most common values are selected from the range from 30% to 70%, depending on the density of the network and node's fluctuation year by year. In general, the selection of parameters should reflect the needs of researchers. For example, one may choose very high value (e.g. 80%) in order to preserve only very similar groups. In another case, it might be necessary to set very low value, e.g. 10% if the network is sparse or the fluctuation is high. In our study, we ran the GED method with alpha and beta parameters varying from 5% to 100%, to see how the number of events varies. Our goal was to have at least one event assigned to each TC. As the splitting and merging events involve several groups, we aimed to have on average slightly more than one event per TC. With this assumption, we selected 30% for both alpha and beta parameters in case of BCN, and 10% for alpha and beta parameters in case of CN. This values produced in total 479 events per 430 groups for BCN, and 492 events per 457 groups for CN.

Correlation between overlap measure and inclusion measure
For BCN, we use the forward and backward intimacy indices to measure the closeness between TCs in consecutive time windows (years). For CN, we considered two types of measure: (i) a simple overlap measure of two groups (the relative fraction of common members), and (ii) an overlap of two groups enriched with the information about the importance of the common members. The latter is suggested by the GED method authors, who named their similarity measure the inclusion measure. One way to evaluate the importance of TC members is to use node centrality measures to rank them within the group. In our work, we are using the Social Position measure 23 (as suggested in the GED method), an idea based on the PageRank algorithm 33 . Saganowski et al. 34 found that using a richer similarity measure allows us to track group evolution more reliably. To better understand the difference between the simple overlap measure and the inclusion measure we compared values obtained with both measures in Figure S1. It turned out that the inclusion measure is on average 20% lower than the simple overlap measure, and the corresponding values, i.e. 30% for the simple overlap and 10% for the inclusion measure, produce roughly the same number of the evolution events. However, the more complex version of the similarity measure (i.e. the inclusion measure), provided slightly better initial prediction results. Therefore, we finally utilized the inclusion measure in our calculations for CN.  Figure S1. The scatter plots for simple overlap measure and inclusion measure for CNs between 1981 to 2010. The left panel (A) is for the alpha parameter, i.e. how the groups in t are close to groups in t + 1. The right panel (B) is for beta parameter, i.e. how the groups in t + 1 are close to groups in t. The sizes of circles are proportional to the number of instances. The red dash lines are y = x for reference only.

18/24
Alluvial diagram for CN Like the bibliographic coupling network (BCN), the co-citation network (CN) can also be visualized in the form of the alluvial diagram. The groups in a CN represent the papers from the past that are coherent and related to a certain topic that stimulates the present research lines. 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 Figure S2. The alluvial diagram of APS papers' references (CNs) from 1981 to 2010. Each block in a column represents a TC extracted from the CN. The height of the block is proportional to the number of papers in the TC. For clarity, only TCs comprising more than 1% of all nodes in CNs are shown. TCs in successive years are connected by streams whose widths at the left and right ends are proportional to the relative overlap percentage. The colours inside a TC represent the relative contributions from different journals.

19/24
The list of features used in the study As we mentioned in Future events prediction, each observation contained 77 features (preselected from the initial 100). The full list of 100 features are showed in Table S1. Many features in this list are proposed for directed social network, therefore are inappropriate for our undirected BCN and CN. The symbol + indicates this feature was used in BCN prediction, while the symbol * indicates this feature was used in CN prediction.

Features group Feature name Feature description
Members/microscopic sum group degree in The sum of indegree 35 of nodes belonging to the community calculated within the community. Indegree is a node measure defining the number of connections directed to the node avg group degree in The average value of indegree of nodes belonging to the community calculated within the community min group degree in The minimum value of indegree of nodes belonging to the community calculated within the community max group degree in The maximum value of indegree of nodes belonging to the community calculated within the community sum group degree out The sum of outdegree 35 of nodes belonging to the community calculated within the community. Outdegree is a node measure determining the number of connections outgoing from the node avg group degree out The average value of outdegree of nodes belonging to the community calculated within the community min group degree out The minimum value of outdegree of nodes belonging to the community calculated within the community max group degree out The maximum value of outdegree of nodes belonging to the community calculated within the community sum group degree total+ * The sum of total degree of nodes belonging to the community calculated within the community. Total degree is the sum of indegree and outdegree avg group degree total+ * The average value of total degree of nodes belonging to the community calculated within the community min group degree total+ * The minimum value of total degree of nodes belonging to the community calculated within the community max group degree total+ * The maximum value of total degree of nodes belonging to the community calculated within the community sum group betweenness+ * The sum of betweenness 35 of nodes belonging to the community calculated within the community. Betweenness is a node measure describing the number of the shortest paths from all nodes to all others that pass through that node avg group betweenness+ * The average value of betweenness of nodes belonging to the community calculated within the community min group betweenness+ * The minimum value of betweenness of nodes belonging to the community calculated within the community max group betweenness+ * The maximum value of betweenness of nodes belonging to the community calculated within the community Continued on next page