GASVeM: A New Machine Learning Methodology for Multi-SNP Analysis of GWAS Data Based on Genetic Algorithms and Support Vector Machines

: Genome-wide association studies (GWAS) are observational studies of a large set of genetic variants in an individual’s sample in order to ﬁnd if any of these variants are linked to a particular trait. In the last two decades, GWAS have contributed to several new discoveries in the ﬁeld of genetics. This research presents a novel methodology to which GWAS can be applied to. It is mainly based on two machine learning methodologies, genetic algorithms and support vector machines. The database employed for the study consisted of information about 370,750 single-nucleotide polymorphisms belonging to 1076 cases of colorectal cancer and 973 controls. Ten pathways with different degrees of relationship with the trait under study were tested. The results obtained showed how the proposed methodology is able to detect relevant pathways for a certain trait: in this case, colorectal cancer.

Introduction
The results of the Human Genome Project [1] and the International HapMap Project [2] made it possible to find genes linked to traits and health problems. Genome-wide association studies (GWAS) have contributed to several new discoveries in human genetics. GWAS exploit the fact that genetic variants that are close together tend to be statistically correlated, something which in genetics is known as linkage disequilibrium [3]. The advances in genome arrays of genetic variations have led to the discovery of many DNA variants associated with complex traits such as those related to diseases.
Nowadays, one of the main criticisms of GWAS is that to date, most of the discoveries have not been applied in a clinical practice [4], but despite this obvious drawback, GWAS do have a great relevance. As an example, it can be said that until the development of GWAS it was not possible to find any gene linked to schizophrenia [5].
Not only are GWAS of interest for the discovery of robust associations, but they also give information about the nature of variations in traits and have contributed to the discovery of new biological knowledge about how DNA variations can affect gene regulation. The variations in the human genome are mainly down to two causes: point

Support Vector Machines
Support vector machines (SVM) are a supervised-learning classification technique that has shown its ability in dealing with classification [16] and regression problems [17,18]. In the case of the present research, SVM is employed for binary classification in cases and controls. Let us suppose they are denoted by {−1, +1}. Therefore, predictors of the form f : R D → {−1, +1} are considered.
The SNPs of each member of the population are represented by a vector x n ∈ R D where D is the number of SNPs employed in the study and the case or control label is given by y n . Given a training dataset consisting of pairs {(x 1 , y 1 ), (x 2 , y 2 ) . . . (x n , y n )} with the objective being to train an SVM model with the lowest classification error. Let x i ∈ R D be a data sample, and consider the function f , R D → R in such a way that x → w, x + b are w ∈ R D and b ∈ R the hyperplane that separates the two classes in the binary classification problem can be written as x ∈ R D : f (x) = 0 .
Please note that w is a normal vector for the hyperplane and b the intercept. When training the classifier, we want to ensure that the examples with positive labels (cases) are on the positive side of the hyperplane, while those with negative labels are on the negative side. These two conditions can be expressed as y n ( w, x n + b) ≥ 0 where y n = −1 or 1.
To find a unique solution, one idea is to choose the separating hyperplane that maximizes the margin between the cases and controls. The margin represents the distance of the hyperplane to the closest examples of cases and controls, respectively, in those cases in which it would be possible to assume that the dataset is linearly separable [19], but this is not the case of the present research.
Let us consider one individual x i , which without loss of generality can be considered as a case and labeled as +1. This case observation should be on the positive side of the hyperplane w, x i + b > 0. The distance of x i to the hyperplane is given by the distance of the orthogonal projection of the point to the plane. Using vector addition, x i can be expressed as follows: In other words, all the case observations must be at least a distance r from the hyperplane. It can be expressed by the following equation: And the optimization problem to be solved can be expressed as: maxr Subject to: y n ( w, x n + b) ≥ r ω = 1, r > 0 Can then be interpreted as the maximization of r while ensuring that the case observations lie on the correct side of the hyperplane. This is usually called the margin maximization parameter.
In these cases, like in the one in the present problem where data are not linearly separable, some examples would fall into the margin region or even on the wrong side of the hyperplane. This model is called soft margin SVM and makes use of the following equations: Subject to: y n ( w, x n + b) = 1 − ξ n , ξ n ≥ 0. ξ n is the slack variable corresponding to each observation. The parameter C > 0 trades off the size of the margin and the total amount of slack that we have, and is called the regularization parameter.
The margin term ω 2 is called the regularizer. The convex duality via Lagrange Multipliers of the previous formula can be expressed as follows: When the Lagrangian is differentiated with respect to the three primal variables, ω, b and ξ, the following result is obtained: Setting the three partial derivatives to zero, it would be possible to express the previous equations by means of their dual negative equations. This converts the problem from a maximization into a minimization one: Although SVM requires linearity for classification, the use of kernels makes this methodology also useful for non-linear problems. A kernel can be defined as a function K : X × X → R for which there is a Hilbert space H and φ : X → H is feature map such that: Kernels must be symmetric and positive semidefinite functions in such a way that the kernel matrix K is symmetric and positive semidefinite. The use of different kernels has different effects on the separating hyperplanes. In this research lineal, polynomial, radial basis functions and sigmoid kernels have all been tested. These kernel functions can be expressed as follows [19]:

Genetic Algorithms
Genetic algorithms can be defined as biologically inspired methods for optimization [20]. The foundations of genetic algorithms can be found in the works of Holland [21], Rechenberg [22] and Schwefel [23].
For their initialization, genetic algorithms require an initial set of candidate solutions for the optimization problem to be solved. Table 1 shows the pseudocode of a genetic algorithm. As can be observed in the table, the first step involves creating an initial population. Data representation and how the initial population is created both have a great importance on the genetic algorithm performance. The second operation performed is the crossover. A non-deterministic crossover function can be defined as X : Ω × Ω → Ω . The result of X x i , x j gives a new population member with the same length as x i and x j and as such, all their elements belong either to x i or x j with a certain probability.
Mutation is employed to inject new strings into the next generation [24], which gives the genetic algorithm the ability to search beyond the confines of the initial population. The mutation function can be expressed as: µ : Ω → Ω . It is like a crossover: a nondeterministic function that assigns to each string member a certain probability of being randomly changed. The fitness computation assigns a value to each member of the population that represents how well they fit to the problem to be solved. Those individuals with the most favorable results of the fitness function are more likely to be selected as parents of the next-generation offspring.
This process is repeated until any termination condition is reached. In the case of the present research, there are certain key parameters to be controlled in order to obtain a fine-tuned version of the algorithm. These parameters are the number of iterations, the population size, crossover and probability of mutation.

The Proposed Algorithm
The algorithm proposed in the present research makes use of both genetic algorithms and support vector machines in order to find out whether a certain pathway, which in this context can be considered in the same way as a set of SNPs, is able to identify cases and controls for a certain trait or illness. Figure 1 shows the flowchart of the proposed algorithm. The first step involves selecting the subset of SNPs that belongs to the pathway under analysis. This means that from the total number of SNPs included in the database, the information required for the analysis is reduced to a selected subset of SNPs of all the members of the population. In other words, the SNPs chosen are only those that belong to the pathway to be studied.
The members of the genetic algorithm (GA) population for this analysis are strings of "1s" and "0s" that indicate which SNPs will form a part of the SVM model to be computed. Please note that "1" means that the SNP will take part of the SVM model and "0" that it will not. In the case of the present research, each member of the GA population has the same length as the number of SNPs that constitute the pathway under analysis. Please note that each GA population has several members, and an SVM model is trained for each one.
All the classification SVM models are trained using as input variables the SNPs with the "1" value and as output, the variable trait that indicates which elements are cases and which are controls. As may be seen in the flowchart, the initial population is formed by rows from an identity matrix selected in a random way up to the completion of the total number of individuals required for the GA population. This means that in the initial population, only one SNP is active in each population member. In other words, it means that after selecting as input information only those SNPs that belong to the pathway under analysis, the initial population is formed by individuals in which only one of those SNPs is active and, afterwards, the different SNPs that belong to the subset that is being employed are switched on and off with the aim of improving the results of the fitness function. The reason for choosing only one SNP in each member of the initial population is that the goal is to get the maximum values of the fitness function while making use of the minimum number of SNPs required and to allow for the importance of every single SNP to be taken into account individually.
In the following populations, the number of SNPs selected in each population member can be more than one as each one of the population members evolves, taking into account genetic algorithms rules in search of the maximization of the value of the fitness function. The fitness function consists of calculating the area under the ROC curve [25] that is obtained when data are classified, making use of the SVM calculated for that member of the population (string of "0s" and "1s") using the active SNPs as independent variables and as dependent variables whether or not the individual suffers from a certain trait, which in the case of the present research is colorectal cancer. In order to avoid problems related to epistasis, members of the population that choose more than one SNP from the same gene have a value of 0 assigned to their fitness function.
When the stop criterion is reached, i.e., the maximum number of cycles allowed, the area under the ROC curve (AUC )value achieved is recorded and the SNPs employed are deleted from the set of SNPs. The process starts again looking for a new SNPs subset for which the AUC value can be as high as possible. The process is then repeated until the SNP set is empty or until a total of 80 cycles have been completed. Afterwards, the same process is repeated 1000 times, making use of permutations of cases and controls labels. Please note that in this field, some of the existing literature recommends 10,000 permutations [26] while in classical literature the number of permutations considered in most of the papers for estimating the power of a permutation test was 1000 [27][28][29] and in some of them only 500 [30,31]. It is also worth noting that previous research [32,33] stated that 1000 is a reasonable number of permutations for a test at the 5% level of significance. Finally, in the field of genomic studies there is software that considers that permutation values from 1000 can feasibly be employed in GWAS [34].  When the stop criterion is reached, i.e., the maximum number of cycles allowed, the area under the ROC curve (AUC )value achieved is recorded and the SNPs employed are deleted from the set of SNPs. The process starts again looking for a new SNPs subset for which the AUC value can be as high as possible. The process is then repeated until the SNP set is empty or until a total of 80 cycles have been completed. Afterwards, the same process is repeated 1000 times, making use of permutations of cases and controls labels. Please note that in this field, some of the existing literature recommends 10,000 permutations [26] while in classical literature the number of permutations considered in most of the papers for estimating the power of a permutation test was 1000 [27][28][29] and in some of them only 500 [30,31]. It is also worth noting that previous research [32,33] stated that 1000 is a reasonable number of permutations for a test at the 5% level of significance. Finally, in the field of genomic studies there is software that considers that permutation values from 1000 can feasibly be employed in GWAS [34].

Design of Experiments
Nowadays, it is well-known that there are no optimal parameter values valid for all problems [20]. This statement is part of the no-free-lunch theorem [35], which states that there is no overall superior optimization algorithm capable of solving every kind of optimization problem.
Parameter tuning strategies treat the parametrization of GA as an optimization problem. In the case of the present research, the GA parameters are tuned with the help of design of experiments methodology (DOE), but there are other possible approaches. For example, the parameter calibration process can be supported on statistical methods [36], a deterministic control whereby an extended scheme is used to control the parameters [37], adaptive control strategies like Rechenberg's mutation rate [22] or self-adaption [20]. In the case of the present research, DOE was chosen.
DOE is a statistical methodology that is employed in order to find relationships among variables that affect certain processes [38]. In other words, DOE makes it possible to see how a simultaneous change in more than one variable will affect the output variable. Although DOE can be applied to both categorical and continuous variables, in the case of the present research all the variables under study with DOE are continuous.
A design of experiments was performed to select the most suitable values of the GA parameters for the algorithm proposed in the present research, using the colorectal cancer pathway as a reference. Furthermore, for designing the experiment, a full factorial design with center points was employed. In a full factorial design, all possible combinations of factor levels are used [39]. For each point, experiments were repeated three times. Please note that each individual experimental setting is referred to as a run and the response measured is called an observation. The present work made use of DOE for the fine-tuning of the GA algorithm. For the DOE analysis, the continuous variables of the research are considered, namely, the number of iterations of the algorithm, the population size and the values of mutation and crossover. The response value is the area under the ROC curve obtained with the SVM model that makes use of the variables selected. Each variable was measured at three different levels and, therefore, a 3 4 full factorial design with 81 experiments, each one repeated three times, was obtained. The variables employed and their corresponding level are shown in Table 2.

Datasets
The database employed in this study belongs to the Colorectal Cancer Transdisciplinary Study (CORECT) project. This was an observational multicentric multi-case control study performed from September 2008 to December 2013. For this research, the subset of information belonging to Leon University Hospital and Hospital of Bellvitge was employed. It contains 1076 cases of, and 973 controls for, colorectal cancer, for which the information from 370,570 SNPs was available.
The cases are incidental and histologically confirmed, with ages between 20 and 85. Those with communication disabilities, physically unable to participate or with a previous diagnosis of colorectal cancer were excluded. For their recruitment, the study staff contacted them at the selected hospitals.
Controls were randomly selected from the population lists assigned to family physicians in the catchment area of the hospitals where the cases were recruited, and with the same sex and age distribution (±5 years). All had been residing in the area of the hospital where the cases were recruited for at least 6 months before and did not present physical or communication impediments.
The protocol of the project was approved by the Ethics Committees of the institutions that took part in the study. The participation of the subjects in the study was voluntary, after signing an informed consent. The confidentiality of the data is guaranteed by eliminating the personal identifiers in the datasets and all the files that include information about the subjects, complying with Spain's Organic Law 15/1999. Please also note that the files have been registered with the Spain Data Protection Agency (Number 2102672171). Access to this information for other researchers is allowed on request.
For the present study, ten different pathways were selected from the KEGG database [40][41][42] so as to include, on the one hand, pathways for which there was already significant scientific evidence that they were associated with the trait analyzed; on another hand, pathways for which the evidence indicated that their association with the trait was improbable, and finally others for which the evidence was inconclusive. Once the ten pathways to be included had been chosen, all the SNPs that belonged to the genes considered in each one of the pathways under analysis were retrieved from the database.

Results
After having fixed the GA parameters to a population size of 5500, 6000 iterations for each cycle and a 100% crossover with a mutation rate of 1%, the proposed algorithm was applied to 10 different pathways.

Design of Experiments
The variables employed in this design of experiments were number of iterations, population size, crossover, and mutation rate. The values tested are presented in Table 2. Each of the combinations of variables combinations was tested three times. Figure 2 shows the main effect plots of these four variables. According to the results obtained, the number of iterations of the GA was fixed at 6000, as there is only a very small increase in the AUC value from 6000 to 8000 (about 0.1%). In the case of the population size, the value of 5500 individuals was considered to be sufficient, as increasing the number of population members to 10,000 only meant an improvement of less than 0.2% in the AUC results. For the mutation rate whose values are presented in logarithmic scale, the maximum of the three values tested was achieved for 1%, which was the center point. Finally, in the case of the crossover rate, 100% was chosen due to it giving the highest performance. For all the pathways analyzed in this research, the number of iterations is 6000, with a population size of 5500, a mutation rate of 1% and a crossover of 100%.

Application of the Algorithm to Different Pathways
After having fixed the GA parameters at a population size of 5500, 6000 iterations for each cycle and a 100% crossover with a mutation rate of 1%, the proposed algorithm was applied to 10 different pathways with different degrees of relationship with the trait under study.
As can be seen in Figure 1, the first step in all the cases involves selecting as part of the initial set only those SNPs of the database taking part in the pathway under study. For

Application of the Algorithm to Different Pathways
After having fixed the GA parameters at a population size of 5500, 6000 iterations for each cycle and a 100% crossover with a mutation rate of 1%, the proposed algorithm was applied to 10 different pathways with different degrees of relationship with the trait under study.
As can be seen in Figure 1, the first step in all the cases involves selecting as part of the initial set only those SNPs of the database taking part in the pathway under study. For example, in the case of the adipocytokine signaling pathway, it has a total of 752 SNPs and in the case of the AMP-activated protein kinase (AMPK) signaling pathway the number is 1812. The initial population of 5500 individuals is always formed by vectors in which only one element is different from zero. It indicates the SNP employed for the SVM classification in the first iteration. Starting from such an SNP, new individuals that combine zeros and ones are formed in search of those SNPs that can provide the maximum AUC value while at the same time taking into account that no more than one SNP from the same gene can be included in the same population member. This iterative process is repeated 6000 times. In the case of the AMPK signaling pathway, the AUC value obtained after this process was 0.584023, and for the adipocytokine signaling pathway it was 0.565382. When the iterations are finished, the SNPs employed in the subset with the maximum AUC obtained are removed and the process is repeated in search of the best remaining SNPs. Although this process would normally be repeated while there were still SNPs pending employment, it was halted after 80 cycles. There are two reasons for this. On the one hand, the time required for these iterations, which is quite high (34.51 s per iteration on average) and on the other hand that although we do not know in advance the exact number of SNPs involved in 80 iterations, in order to compare the results of pathways of different lengths, a number of repetitions was chosen that would be feasible for any of the pathways.
After having run the algorithm for all the pathways included in the study, the same process was repeated 1000 times for each one but permuting phenotypes of cases and controls while preserving the total number of 1076 cases of colorectal cancer and 973 controls. The results obtained were compared.
The main results are detailed in Table 3. This table presents the total number of SNPs that are included in each of the 10 pathways under analysis. As was mentioned before, the algorithm was repeated 80 times in all cases. This means that not all the SNPs were employed for the process of classification of cases and controls. For example, for the adipocytokine signaling pathway, which has a total of 752 SNPs, only 496 were employed in any of the iterations for the classification of individuals in cases and controls. Figure 3 shows the AUC values of the 80 iterations performed for the adipocytokine signaling pathway in the case of cases and controls (phenotype) and for 5 different permutations of the 1000. For the graphical representation and for the comparison of wins in pathway versus permutated, AUC values are ordered from highest to lowest. In this case, the phenotype curve (in green) does not seem to classify cases and controls in a better way than the permutated ones. Please note that as can be observed in Table 3, the average AUC value of the 80 cycles of phenotypes is 0.535858, while in the case of case of permutated cases and controls it is 0.537543, which means it is 0.31% lower. The column called win subsets indicates the percentage of times when the AUC value of the phenotype is higher than the permutated ones.
Something similar happens in the case of the insulin resistance pathway, where the percentage of the win pathways is 29.75% and the values of the AUC and the average permutated AUC are 0.555483 and 0.556201, respectively (−0.13%). Therefore, in this case, as in the adipocytokine signaling pathway, there does not seem to be any relationship between the two pathways and colorectal cancer. Please also note how in Figure 4 the curve of Phenotype does not seem to be higher than the permutated ones. In addition, for the longevity regulating pathway, whose curve is represented in Figure 5, the situation is similar and no significant influence of those pathways on colorectal cancer can be reported. Table 3. Pathways under analysis. Total number of single nucleotide polymorphism (SNPs) per pathway (Tot. SNPs), SNPs employed in the 80 iterations by the non-permutated phenotypes (SNPs employed), average area under the receiver operation curve (AUC) obtained in the 80 iterations by the non-permutated phenotypes (AUC), AUC obtained by the permutated phenotypes (AUC perm.), percentage of non-permutated AUC values than are higher than the maximum permutated AUC value (win subsets).  Something similar happens in the case of the insulin resistance pathway, where the percentage of the win pathways is 29.75% and the values of the AUC and the average permutated AUC are 0.555483 and 0.556201, respectively (−0.13%). Therefore, in this case, as in the adipocytokine signaling pathway, there does not seem to be any relationship between the two pathways and colorectal cancer. Please also note how in Figure 4 the curve of Phenotype does not seem to be higher than the permutated ones. In addition, for the longevity regulating pathway, whose curve is represented in Figure 5, the situation is similar and no significant influence of those pathways on colorectal cancer can be reported. Something similar happens in the case of the insulin resistance pathway, where the percentage of the win pathways is 29.75% and the values of the AUC and the average permutated AUC are 0.555483 and 0.556201, respectively (−0.13%). Therefore, in this case, as in the adipocytokine signaling pathway, there does not seem to be any relationship between the two pathways and colorectal cancer. Please also note how in Figure 4 the curve of Phenotype does not seem to be higher than the permutated ones. In addition, for the longevity regulating pathway, whose curve is represented in Figure 5, the situation is similar and no significant influence of those pathways on colorectal cancer can be reported.   The opposite case applies for apelin signaling, mitochondrial biogenesis and colorectal cancer. In these three cases, the average AUC value obtained for cases and controls are clearly higher than for the permutated phenotypes. In the case of the apelin signaling pathway, it is 5.15% higher in the case of phenotype when compared with the permutated solutions. In the case of mitochondrial biogenesis, it is 3.23% and for the colorectal cancer pathway, the value is 2.45%. In all the cases, the AUC values obtained are higher in the phenotypes than in the permutated cases with 100% of winning cases for apelin signaling pathway, colorectal cancer pathway and mitochondrial biogenesis. Figures 6-8 clearly show how the phenotype curves are higher than the permutated ones. Although in the case of Figure 9, where the AMPK signaling pathway is represented, it does not seem to be as clear as in the three previous cases, the AUC value is 2.26% higher when compared with permutated cases. Please note that in 89.75% of cases the values obtained are higher in the phenotypes than in the permutated cases which, from our point of view, would mean that there is certain influence of this pathway on the trait under analysis.  The opposite case applies for apelin signaling, mitochondrial biogenesis and colorectal cancer. In these three cases, the average AUC value obtained for cases and controls are clearly higher than for the permutated phenotypes. In the case of the apelin signaling pathway, it is 5.15% higher in the case of phenotype when compared with the permutated solutions. In the case of mitochondrial biogenesis, it is 3.23% and for the colorectal cancer pathway, the value is 2.45%. In all the cases, the AUC values obtained are higher in the phenotypes than in the permutated cases with 100% of winning cases for apelin signaling pathway, colorectal cancer pathway and mitochondrial biogenesis. Figures 6-8 clearly show how the phenotype curves are higher than the permutated ones. Although in the case of Figure 9, where the AMPK signaling pathway is represented, it does not seem to be as clear as in the three previous cases, the AUC value is 2.26% higher when compared with permutated cases. Please note that in 89.75% of cases the values obtained are higher in the phenotypes than in the permutated cases which, from our point of view, would mean that there is certain influence of this pathway on the trait under analysis. The opposite case applies for apelin signaling, mitochondrial biogenesis and colorectal cancer. In these three cases, the average AUC value obtained for cases and controls are clearly higher than for the permutated phenotypes. In the case of the apelin signaling pathway, it is 5.15% higher in the case of phenotype when compared with the permutated solutions. In the case of mitochondrial biogenesis, it is 3.23% and for the colorectal cancer pathway, the value is 2.45%. In all the cases, the AUC values obtained are higher in the phenotypes than in the permutated cases with 100% of winning cases for apelin signaling pathway, colorectal cancer pathway and mitochondrial biogenesis. Figures 6-8 clearly show how the phenotype curves are higher than the permutated ones. Although in the case of Figure 9, where the AMPK signaling pathway is represented, it does not seem to be as clear as in the three previous cases, the AUC value is 2.26% higher when compared with permutated cases. Please note that in 89.75% of cases the values obtained are higher in the phenotypes than in the permutated cases which, from our point of view, would mean that there is certain influence of this pathway on the trait under analysis.   The last three pathways under study were glucagon signaling (Figure 10), Huntington's disease ( Figure 11) and insulin signaling ( Figure 12). In these three pathways, in a similar way to the case of the AMPK signaling pathway, the AUC value of the phenotype is slightly higher than the average value obtained for the permutated ones. Additionally, in most of the cases (from 82.50% to 96.50%), the AUC obtained in the phenotype iteration is higher than in the permutated one.  The last three pathways under study were glucagon signaling (Figure 10), Huntington's disease ( Figure 11) and insulin signaling ( Figure 12). In these three pathways, in a similar way to the case of the AMPK signaling pathway, the AUC value of the phenotype is slightly higher than the average value obtained for the permutated ones. Additionally, in most of the cases (from 82.50% to 96.50%), the AUC obtained in the phenotype iteration is higher than in the permutated one.  The last three pathways under study were glucagon signaling (Figure 10), Huntington's disease ( Figure 11) and insulin signaling ( Figure 12). In these three pathways, in a similar way to the case of the AMPK signaling pathway, the AUC value of the phenotype is slightly higher than the average value obtained for the permutated ones. Additionally, in most of the cases (from 82.50% to 96.50%), the AUC obtained in the phenotype iteration is higher than in the permutated one.

Pathway
In summary, taking into account the results obtained with the algorithm proposed in the present research, it can be said that there is a clear relationship linking apelin signaling, colorectal cancer and mitochondrial biogenesis pathways with colorectal cancer. A weak relationship with colorectal cancer was found for AMPK signaling, glucagon signaling, Huntington's disease and insulin signaling pathways. Finally, no relationship with colorectal cancer was found for adipocytokine signaling, insulin resistance or longevityregulating pathways.  In summary, taking into account the results obtained with the algorithm proposed in the present research, it can be said that there is a clear relationship linking apelin signaling, colorectal cancer and mitochondrial biogenesis pathways with colorectal cancer. A weak relationship with colorectal cancer was found for AMPK signaling, glucagon signaling, Huntington's disease and insulin signaling pathways. Finally, no relationship with colorectal cancer was found for adipocytokine signaling, insulin resistance or longevityregulating pathways.

Discussion
From the authors' point of view, the results obtained are coherent with previous results to be found in existing literature. Regarding the AUC values, it must be stated that although they would seem to be low, they are in line with other research [43]. In the case of the AMPK signaling pathway, previous studies found that AMPK promotes the survival of colorectal cancer stem cells [44]. This study found that colorectal cancer stem cells show a higher level of antioxidant genes and have a lower level of reactive oxygen species than non-colorectal cancer stem cells. This would be because the colorectal cancer stem cells also possess more mitochondria mass and show higher mitochondrial activity. In the case of this study, higher AMP-activated protein kinase (AMPK) activity was observed in these colorectal cancer stem cells.
Another study that included patients with stage II or III of colorectal cancer, amongst whom the 5-year survival rate was between 50% and 87%, found that the AMPK encoded

Discussion
From the authors' point of view, the results obtained are coherent with previous results to be found in existing literature. Regarding the AUC values, it must be stated that although they would seem to be low, they are in line with other research [43]. In the case of the AMPK signaling pathway, previous studies found that AMPK promotes the survival of colorectal cancer stem cells [44]. This study found that colorectal cancer stem cells show a higher level of antioxidant genes and have a lower level of reactive oxygen species than non-colorectal cancer stem cells. This would be because the colorectal cancer stem cells also possess more mitochondria mass and show higher mitochondrial activity. In the case of this study, higher AMP-activated protein kinase (AMPK) activity was observed in these colorectal cancer stem cells.
Another study that included patients with stage II or III of colorectal cancer, amongst whom the 5-year survival rate was between 50% and 87%, found that the AMPK encoded in gene α1 was overexpressed in patients who suffered from colorectal cancer. For its authors, the AMPK encoded in gene α1 regulate the glutathione reductase (GSR) phosphorylation, possibly through residue Thr507, which enhances its activity. They suggested the suppression of AMPK expression in gene α1 by using nano-sized polymeric vector to induce a favorable therapeutic effect.
In the case of the apelin signaling pathway, which was also found to be relevant to colorectal cancer, there are some studies linking it to colorectal cancer [45]. Apelin is an endogenous ligand of the apelin receptor (APJ), a seven-transmembrane G protein-coupled receptor [45]. It can be found in the brain and also in peripheral organs like the heart, the lungs, blood vessels, and adipose tissue. It is involved in regulating cardiac and vascular function, heart development, and vascular smooth muscle cell proliferation. According to previous research, apelin is not only related to colorectal cancer, but also to others like lung cancer, gastroesophageal, hepatocellular carcinoma, prostate cancer, endometrial cancer, oral squamous cell carcinoma, brain cancer, and tumor neoangiogenesis. This means that Apelin/APJ may be a potential anticancer therapeutic target. A study suggested that the APJ receptor antagonist F13A significantly reduced cellular proliferation [46]. Another study [47] found that Apelin receptor is co-expressed in colorectal cancer cell lines and its activation leads to adenylyl cyclase inhibition and Akt phosphorylation. For the authors of that research, apelin and its receptor might be co-expressed in the tumor compartment where this co-expression would underlie a constitutive activation of apelin signaling and create a functional autocrine loop. It was the first study that reported that apelin peptide is highly expressed in human colon adenomas and tumors [47]. This co-expression was also observed in several colorectal cancer cell lines. In the LoVo cell line, quantitative real-time polymerase chain reaction (qRT-PCR) experiments and apelin-induced Akt phosphorylation confirmed the concomitant expression of both ligand and receptor. In addition, apelin behaved as an anti-apoptotic peptide, by reversing caspase activation and poly ADP ribose polymerase protein (PARP) degradation induced by theMG132 proteasome inhibitor. Another study [48] that measured apelin, and its receptor mRNA, and protein expression levels in tumor tissue of 56 surgically treated colorectal adenocarcinoma patients and compared them with 27 healthy controls, found that serum levels of apelin and its receptor were increased in colorectal cancer patients in comparison to controls, which leads to the conclusion that apelin could be an important factor in the progression of colorectal carcinoma. The finding of the colorectal cancer pathway as being significant by our algorithm is not surprising, as it can be considered as the pathway of reference.
Mitochondria are semiautonomous organelles that participate in energy metabolism, free radical production, and apoptosis. Apart from the nucleus, the mitochondrion is the only cellular organelle that contains its own genome and genetic machinery [49,50]. Mitochondrial biogenesis is an essential process by which new mitochondria are obtained, and is one which requires coordination between the nuclear and mitochondrial genomes [51]. Mitochondria, as well as most of the processes related to them, are closely linked to the genesis of cancer [52]. For this reason, it is essential to study mitochondrial biogenesis, as well as to find out what happens to these organelles during tumor processes. The progression of CRC in humans is closely linked to mitochondrial alteration, increased production of free mitochondrial oxide radicals, and oxidative stress [53].
An article in a literature review article [54], aimed at evaluating whether increased or decreased peroxisome proliferator-activated receptor gamma coactivator 1-α (PPARGC1A or PGC1α) expression affects the development of colorectal cancer, found that an altered expression of PGC1α modifies colorectal cancer risk and mitochondrial biogenesis is regulated by PGC1α. According to this study, it seems plausible that the proposed algorithm found a relationship between the mitochondrial biogenesis pathway and colorectal cancer.
Glucagon increases the production of glucose by increasing glycogenolysis and gluconeogenesis in the liver, and by reducing glycogenesis and glycolysis. The release of glucagon in response to food consumption depends on the type of meal that has been eaten. If a meal is rich in carbohydrates, blood glucagon levels fall to prevent an undue rise in the level of circulating glucose. Conversely, when a protein-rich meal is eaten, the blood glucagon level rises. Nowadays cancer is known to be one of the major causes of death in diabetic patients, and an association between antidiabetic drugs and the risk of cancer has been reported [55]. Glucagon is nowadays recognized as a pivotal factor implicated in the pathophysiology of diabetes. A recent study has found [56] expression of the glucagon receptor in colon cancer cell lines and in colon cancer tissue obtained from patients. According to this study, glucagon significantly promoted colon cancer cell growth. Molecular assays showed that glucagon acted as an activator of cancer cell growth through deactivation of AMPK and activation of mitogen-activated protein kinase (MAPK). Another study [57] found the relationship between glucagon signaling pathway and endometrial cancer.
A study published in 2002 [58] stated that Huntington's disease provides clues about cancer and that it would be a marker of certain cancers like colorectal cancer. It can be said that, in general, people with Huntington's disease have been observed to have lower rates of cancers [59] and although the relationship of Huntington's disease with prostate cancer has been reported [60], no similar study has been found for colorectal cancer.
The insulin signaling pathway is another of the pathways where a moderate relationship with colorectal cancer was found. It has been reported that the modification in the individual values of plasma insulin levels due to diet may affect the risk of suffering from colorectal cancer [61]. A similar result was found by other researchers in a study performed with a sample of postmenopausal women [62]. Although there are many studies in this line [63,64] and it is known that genetic variants in metabolic signaling pathways may interact with lifestyle factors such as dietary fatty acids influencing colorectal cancer risk, these interrelated pathways are not fully understood [65].
As mentioned before, the proposed algorithm did not find any relationship linking adypocitokine signaling, insulin resistance and longevity-regulating pathways with colorectal cancer. In the case of the adipocytokine signaling pathway, no relationship has been found in the existing literature with colorectal cancer. However, its relationship with atherosclerosis, diabetes and breast cancer [66] has been reported. In the case of insulin resistance, as far as the authors know, the relationship is not clear enough at this point in time [67]. Finally, in the case of the longevity-regulating pathway, there were no studies found linking it with colorectal cancer.

Conclusions
This paper presents a novel method called GASVeM, which is based on two wellknown machine learning methodologies-genetic algorithms, and support vector machines. Although the results achieved appear promising, as usually happens in machine learning methodologies applied to GWAS, it is difficult to find a direct biological link between the SNPs involved in the results and the trait under study. In spite of this, through studying existing literature it has been possible for the authors to find previous well-known relationships between the relevant pathways and the trait under analysis.
From the authors' point of view and due to a lack of a machine learning gold standard for GWAS analysis, the present method could be of interest for future GWAS. In this direction, based on the results obtained, it would be interesting, on the one hand, to work on studying the classification capacity when information from several pathways is combined and, on the other hand, to try to replicate the results obtained in other databases and in the analysis of other pathologies, to validate the usefulness of the method under different conditions of use.
We also agree with those authors that consider that we are in the infancy of the use of machine learning in GWAS [68] as we are still quite a long way from achieving gold standard methods producing consistently validated biological insights. In addition, we would like to highlight the characteristics of this kind of database, whereby a high number of SNPs (columns) when compared with the number of cases (rows) cause a kind of problem with GWAS that is difficult to deal with from a machine learning point of view and that, in the case of the present research, is present in the need for an a priori SNPs selection based on pathways. It is the authors' opinion that this problem has a great impact on the reproducibility of results when the same algorithm is applied to a different database.
Finally, it must be said that we are aware that the translation of the results obtained with this method to a population-based clinical practice to carry out a personalization of interventions based on genomic data still requires further steps to be taken before the selection capacity can be refined. However, we consider that the method has demonstrated, as was our objective, a good ability to discriminate which pathways are associated with the event and which are not, through the choice of a limited set of SNPs. We consider the ability to classify individuals to be a second step to be taken in the development of the model, through lines of study such as the inclusion of several pathways in the model. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.


Mistake Captioning: A Machine Learning Approach For Detecting Mistakes and Generating Instructive Feedback

Giving feedback to students is not just about marking their answers as correct or incorrect, but also finding mistakes in their thought process that led them to that incorrect answer. In this paper, we introduce a machine learning technique for mistake captioning, a task that attempts to identify mistakes and provide feedback meant to help learners correct these mistakes. We do this by training a sequence-tosequence network to generate this feedback based on domain experts. To evaluate this system, we explore how it can be used on a Linguistics assignment studying Grimm’s Law. We show that our approach generates feedback that outperforms a baseline on a set of automated NLP metrics. In addition, we perform a series of case studies in which we examine successful and unsuccessful system outputs.

Introduction
Giving feedback is one of the most critical parts of education and training. It allows the instructor to correct errors in students' understanding and guide them to the correct solution. In automated learning systems, such as intelligent tutoring systems (ITSs) the feedback and hints are focused on getting the student to complete the assignment, but this may not take into account their mastery of that skill. While better feedback can be made by considering student intent, we find that current models are lacking in this regard. They are typically focused on learning objective mastery and may not necessarily give the whole picture when determining why students answered something the way they did.
But there is a natural solution to this, in that skilled educators are capable of intuitively modeling a student's thought process and determining mistakes in it. This skill is often developed over time. The difficulty of this task for educators depends on the field of study, with some more free form fields like writing and programming being particularly difficult. By indirectly modeling a student's intent by instead modeling how the instructor gives feedback, we can use this model to generate novel feedback sentences. We propose a system that when given a mistake that has been identified in the task of deriving sound changes from Proto-Indo-European to Proto-Germanic, will automatically generate an explanation that identifies the mistake and the reasoning to why this sound change should not apply in this case. We call this method Mistake Captioning as it takes inspiration from image captioning methods and post-hoc AI rationale generation methods such as Ehsan et al. (2019) Automated Rationale Generation.
We apply a suite of automated machine translation and image captioning metrics to test our system. We also take a closer look at 3 cases from a more complete retraining of the network to better understand how the network is performing and how well the metrics perform regarding our task. We find that this method is a good start, but more data and work need to be done to be able to integrate this into a serious learning environment.

Background and Related Work
In the space of ITSs a problem is how to get a student to the correct solution. One of the ways to achieve this is by having an expert to author the correct solution and the path to take to get there (Marwan et al., 2020;Unnam et al., 2019;Ariely et al., 2020). These all require the use of an expert to either structure the assignment such that feedback can be extracted or to label assignments to create a system that can learn the expert knowledge. In some domains this is not necessary as the students can provide the data for the path to the solution themselves. This data-driven approach uses data from multiple students to create all possible paths that lead to a solution, as is shown by Stamper et al. (2008) in their work on Hint Factory. This was later extended to create paths that optimize for productivity metrics by Maniktala and Barnes (2020). In vast solution spaces generating all possible paths to the solution is not feasible, and in Rivers and Koedinger (2017) work this is done by abstracting states so that a single state can represent a large space of states.
Our method does not attempt to establish a correct answer like these previous methods as they do not tackle an underlying problem as to why students make mistakes. We structure our problem not as a method of guiding the student to the correct path, but by trying to correct the thought process that led them to their current answer. This is similar to the methods used in image captioning (You et al., 2016) where they train a model to recognize what is in an image by relating human authored captions to images. This model is then used to generate captions on similar features in new pictures. This is also like the work done in post-hoc AI rationale generation like that of Ehsan et al. (2019). In this work they relate a set of human actions to human rationales of those actions for the use of explaining AI behavior. In our system, we make use of human actions in the form of mistakes and relate them to human explanations of the mistake.

The Proto-Indo-European Language and
Grimm's Law Our application area for this paper is on an assignment concerning Grimm's Law, which is a series of sound changes that occurred in the evolution of Proto-Indo-European (PIE) into Proto-Germanic (PGmc). PIE is a reconstructed language that attempts to recreate the common ancestor of Indo-European language family. While no direct evidence remains of PIE, the similarities between the languages in this family indicate that they come from a common ancestor. Despite these similarities, there are still differences in the languages that can be attributed to shifts in pronunciation as time went on and communities speaking this language became isolated from others. PGmc. is also a reconstructed language, serving as the source of all Germanic languages. Sound changes can be traced from PIE to PGmc. We focus our task on a set of sound laws that describe how stop consonants change from PIE to PGmc collectively known as Grimm's Law (Campbell, 2013) and also include instances where Grimm's Law does not occur. Grimm's Law consists of three different changes that all contribute towards a single shift known as a chain. These changes are as follows: 1. PIE voiceless stops change to voiceless fricatives.
2. PIE voiced stops change to voiceless stops.
3. PIE voiced aspirated stops change to voiced stops or fricatives.
Together they are chained such that affected sounds in PIE are shifted one step to their form in PGmc. This limited set of shifts is what is represented in our data and is the reason for some of the limitations.

Data
The data used in this paper is centered around the linguistics PIE rule, Grimm's Law (Campbell, 2013). Since the task is to generate a reasoning for a mistake made by students, we have opted to organize the data as a fill-in-the-blank task, as one would show up as an assignment. Since we are focusing on captioning the mistake, the changes from PIE to PGmc must contain a mistake and an explanation for that mistake. To create these erroneous responses, we had a linguistics expert generate a number of entries fulfilling these categories: a PIE word; a PGmc form of that word with blanked areas demarked by underscores; a faulty PGmc word that contains a mistake; and an explanation for what that mistake was and how to fix it. The PIE word is the original form of the word that the student is tasked to change. The blanked PGmc form is blanked only in the areas that the student would be tasked to change. There may be multiple blanks if there are multiple spots where the change occurs. The faulty PGmc form contains only a single mistake each, so that the explanations can be more focused on that specific mistake. When multiple mistakes are identified for each PIE word, they are separate entries with separate explanations. Likewise, if there are multiple blanks to fill in for a word, then each blank will have a separate set of entries for its possible mistakes. The explanations were written with a simplified style as to not introduce too much detail into the explanations that

PIE Blanked PGmc
Faulty PGmc Explanation pisḱós is az fishazḱ ultimately does become h, but not when immediately following stops & fricatives pisḱós is az biskaz p can shift to b, but only when the middle of the word preceding an accented vowel pisḱós is az fisgazḱ can shift to g, but only when preceding an accented vowel in the middle of a word pisḱós is az fisxazḱ normally changes to x, but not when preceded by a stop or fricative Table 1: The data entries for the mistaken sound changes in the PIE word for "fish" resulted in each explanation being a single sentence following a general format. This sentence is in the rough form of "the mistake that was made, reason why it is not applicable in this case". For example, in Table 1 the explanation for why fishaz is wrong hasḱ ultimately does become h as the mistake that was made and but not when immediately following stops & fricatives as the reason why it is not applicable in this case.
This data was gathered by a single Linguistics Expert, with minor input from the authors as to style and structure to better suit this data for machine learning. Our linguistics expert is an associate professor of linguistics that frequently teaches classes on PIE. There were no attempts to control for consistency in tone, tense, or voice style, so such variations do occur in the data. Because the sentences follow a general format, there are some exact explanation matches in the data. Since this was the application of a small set of rules over a larger set of words, finding common letter changes and explanations could not be avoided. We collected 163 entries of PIE words, fill-in-the-blanks, faulty PGmc words, and their corresponding explanations. These came from 55 unique PIE words and contained 53 unique explanations to cover all cases found in the data. The most common explanation occurred 15 times, while others only occurred once in the entire dataset. Even those that occurred once often shared similar phrasing with other explanations since they were drawn from a common set of rules and a common set of consonants.
The data was further structured for training as seen in Table 2. To preserve the question that was being hypothetically asked to solve and the answer, we combined the PIE form and the faulty PGmc form, separated by a space, as our model input. For our model output we used the explanations without adding anything, but we did remove and change some characters. We removed all instances of ending punctuation like periods and exclamation points and spaced out commas to count them as separate tokens. The white space was also normalized and any leading or trailing white space was removed. While not shown in the table, the input was tokenized character by character to help preserve the differences between the faulty PGmc, while the output was done word by word.

Data Noise
Our data is very limited in scope due to the effort it requires to generate it. If we were to train our model only using the data as it is, we would likely run into issues with overfitting and overtraining. To alleviate this, we have injected noise into the data to be able to train longer. For our training output, creating noise is a straightforward process. In every training iteration, the output has 30% of its words masked out. This value was chosen to keep most of the explanation intact while still having a significant amount noised as overfitting is a serious concern. This means that at no training iteration has the model ever seen a complete output sentence, but with enough training iterations the unmasked portions should have overlapped enough to reveal the complete sentence.
The input was a more complicated process. Since the input is smaller than the output, we made the attempt to preserve some of the important parts, namely the actual changes from PIE to PGmc. Since these are fill-in-the-blank tasks, we have access to the parts that contain no mistakes and using this part of the dataset we identified the portions

Input
Output Maskable Indices pisḱós fishazḱ ultimately does become h , but not when immediately following stops & fricatives [10,11,13,14] pisḱós biskaz p can shift to b , but only when the middle of the word preceding an accented vowel [10,11,13,14] pisḱós fisgazḱ can shift to g , but only when preceding an accented vowel in the middle of a word [10,11,13,14] pisḱós fisxaz k normally changes to x , but not when preceded by a stop or fricative [10,11,13,14]  of the Germanic words that are considered safe to mask out in the form of maskable indices. These are shown in the 3rd column of Table 2, though these may not visually line up to the index of the words since letters that are accented are encoded as multiple characters. Since this only affected the PGmc part of the already smaller input, instead of masking 30% of the characters we only mask a single character at a time and limit it to only 30% of the time. This hopefully allows us to make a more robust set of training examples without disrupting the underlying meaning too much.

Model
For our model we use a sequence-to-sequence network. Sequence to sequence networks utilize two recurrent neural networks: an encoder and a decoder. The encoder encodes the sequential data into a fixed length context vector that is meant to represent the important elements of the input. This context vector is then used by the decoder to generate a natural language output. In our case, this output is an explanation, but this can be used for other tasks. Sequence to sequence networks have been used for translation (Sutskever et al., 2014) and post-hoc rationale generation (Ehsan et al., 2019) tasks, and while our task is not a traditional translation or rationale generation task it can be modeled as one. The network learns to take the input sequence of our PIE and faulty PGmc words as a character sequence and translate them into a word sequence for our explanations. By connecting the input to the explanation, we are learning humanlike reasoning for why the mistakes happened. We use Gated Recurrent Units (GRUs) to take advantage of their faster training time and better use with smaller data (Chung et al., 2014). We also make use of an attention mechanism (Graves, 2013) on the decoder to learn to focus on the important parts of the output explanation. The model was trained on 75,000 iterations using the same 163 question and explanation pairs, but as mentioned previously in each iteration some characters and words were masked to prevent the model from seeing them. During each training iteration a random sample from the data is run through the network after having noise applied.

Experiments
We conducted two experiments to evaluate our method: A quantitative approach using a set of automated metrics to judge the candidate explanations, and a qualitative case study to demonstrate what types of outputs were generated. For both the qualitative and quantitative evaluation we used 5 different automated metrics: BLEU (Papineni et al., 2002), METEOR (Banerjee and Lavie, 2005), ROUGE (Lin, 2004), CIDEr (Vedantam et al., 2015), and BERTScore (Zhang et al., 2019).
For our quantitative evaluation, we used k-fold cross-validation to control for the variance in our data. We compare our approach against a majority baseline. Since several of the explanations in the data were repeated, we decided that a better baseline would be the majority baseline and used the most repeated explanation in our data. This explanation was used 15 times in the data, and for the majority baseline we set the generated output to all be this exact sentence. This explanation can be seen in Table 4 in the baseline incorrect, candidate sentence. While the majority baseline has a low accuracy, about 9.2%, it was chosen to show that the method is viable before continuing to a human evaluation.
For our qualitative analysis, we opted to look at outputs for a retrained model that was trained on all but 8 random examples of the data. Since our dataset is small, we decided that training on all but a small number of examples would give better insights on what the network was learning. We examined how the model performed in generating explanations for each of the 8 examples that were held out during training. From these 8, we used their automated metric performance in conjunction with their correctness to divide them into three categories: correct, partially incorrect, and incorrect. We then explore why this occurred and reason about the implications of these results in practice.
During training, we use the following hyperparameter values for our model. We used a teacher forcing ratio of 0.5, and a learning rate of 0.01. For the hidden layers (consisting of the attention and GRU layers), a fixed size of 256 was used. All outputs were limited to a maximum length of 30 but this limit was not reached. To test the model, we implemented beam search to select the best candidate explanation instead of the default greedy search. In addition to improving the accuracy of the model this allowed us to look at multiple final outputs, each scored by the model. This is useful in manually judging if the alternatives to the top output were more correct. We limited the beam width to 5.

Results
We found that our method consistently outperformed the baseline in all metrics as can be seen in Table 3. This provides evidence that our model is working to learn when to apply and how to recreate humanlike feedback in response to mistakes. We also found that the leave-8 retraining performed better than the 5-fold models. This trend can be visually seen in Figure 1.
Applying a 2-sided t-test we find that the difference between the baseline and the 5-fold is statistically significant (p < 0.0001) in all cases. When comparing the baseline to 8 case studies this also holds true (p < 0.001). There is also no statistically significant difference between the case studies and the 5-fold.

Discussion
In this section, we discuss the results in more detail, going over the results of the automated metrics and case studies.

Automated Metrics
While the automated metrics do show that our method outperforms our baseline, they are not perfect. One finding that came out of this evaluation is that automated metrics can be misleading when evaluating performance. This is because automated metrics for evaluating language generation systems often measure how well generated sentences overlap with reference sentences. The change of a single word alters the semantic meaning of an explanation, however, which can fool the metrics into scoring it higher than it is. Likewise, a generated explanation can be nearly correct but contain little similarity in structure to the original. This is ultimately a limitation of the data since we only have a single reference explanation for each faulty PGmc example. In the data there exist explanations that explain a concept similar to other explanations but are differently structured. If the model lifts the structure of one explanation to the correct change and reasoning of a change, there is no guarantee that this will match the reference explanation and will be scored more poorly. This does not mean that these metrics are not useful, though, just that the limited scope of our data makes it less suited to these metrics. There is a general trend between how well the explanation performs compared to how correct it is, as we will see in the case studies. In general, this seems to suggest that how  well these metrics will perform will change if we get more data, both on the number of PIE words and the number of possible explanations for each mistake.

Case Studies
To take a better look at the generated explanations we have a separate leave-8-out training. In these 8 that we left out and tested on we found three categories: correct, partially incorrect, and completely incorrect. Correct explanations include the correct mistaken change and the correct reason it was not applicable, while partially incorrect only include one of the two. Incorrect explanations contain neither. Of the 8, 4 of these where correct, 3 were partially incorrect, and only one was completely incorrect. We have selected 3 examples to match these 3 categories to look more closely at.

Completely Incorrect
There was only a single example of the completely incorrect case, which is shown in Incorrect case in Table 4. This explanation does not identify the sound change givingḱ instead of t for the sound change. For the reasoning it seems like it is on the correct track by mentioning preceded by stop or fricative but reverses whether the change occurs because of this. In fact, this is a case of the candidate explanation having no relation to the reference explanation, nor is close to any reference explanation for that PIE word, nor the PIE and PGmc input. The candidate explanation refers to aḱ which does not appear in input but, as will be expanded on later, this may be a case of encoding. Due to the way that the input is interpreted on a character-bycharacter basis the model may have recognized the presence of the k and conflated it with aḱ. It also used the explanation that is most repeated in the dataset, which suggests that it picked the most popular explanation as opposed to correctly learning this change.

Partially Incorrect
For the partially incorrect case we find that the model could identify something correct but not everything. Sometimes this was only identifying the correct sound change or the correct reasoning, but in the case that we are focusing it generated a completely new explanation that was one word from being correct. This candidate explanation changed the word must to only which completely changes the meaning of the explanation. In our explanation generation we saved the top 5 sentences, and the correct explanation does appear as the second. This may be another encoding issue since the phrase k only appears several times in our dataset. We used beam search to attempt to avoid this problem, but in this case, it does not seem to have been enough. The complete phrase that was generated does not appear as a reference explanation in any of the examples, so it is unknown as to why it chose this one. It is likely that having more reference explanations would help as it would force the model to learn more important features.
In another case of partially incorrect, we found that the exact reference explanation that never appeared in that training data, which means that the model could not have possibly generated a completely matching explanation.

Correct
The correct category is marked by when the model was correctly able to generate an explanation for the mistake that was made, matching both the change and reasoning. This is also the largest category, comprising 4 of the 8 examples. These score highly on all the metrics. This category has sound changes that are well represented in the training data, with no reference explanation appearing less than 5 times. In the case shown in Table 4  reference explanation appeared 8 times which is relatively high, but not unique to this category. The completely incorrect case had its reference sentence appear 12 times in the data, meaning that this is likely not the most prominent factor in whether a generated explanation will be correct or incorrect. In all these cases we find that the automated metrics to judge these explanations are generally higher when the answer is more correct, and generally lower when it is more incorrect. This holds in cases when there is only partial correctness. As mentioned before this is likely due to our limited data and may break if we have more reference sentences. We hope to expand the data to cover other changes in linguistics, though the method should be able to be transferring to other domains too.

Future Work
Due to the way that our method works in identifying a single mistake at a time, we can use the fill-in-the-blank task to generate explanations on multiple mistakes. This is done by first comparing the answer given to the correct answer and finding the difference between the two and isolating each individual mistake. Each of these individual mistakes can be then applied to the original correct answer to create separate faulty answers which can subsequently be run through our system. Using this method, we can also change the fill-in-the-blank task to a simple response, where the student is tasked to correctly make all the changes to the PIE word to produce a PGmc word.
Either through separate models for separate rules (and enough annotations on the input) or by training a single more complex model, it may be possible to create a system that can caption all mistakes for a given task. This could work in conjunction with automatically grading the assignments to quickly provide feedback on the mistakes.

Conclusion
In this paper we show a method of the novel task of creating automated captions of mistakes. These captions serve as explanations for what the thought process behind the mistake was, and why it is not applicable in this case. We apply two methods to test our model, comparing against a baseline with a suite of automated metrics, and manually identifying and analyzing the sentences in a series of case studies. Our experiments show that our method has promise in creating these automated captions but that there are significant challenges that need to be overcome.
We set out tasks in the domain of Linguistics pedagogy, and in this we can make improvements with representation of the input data and with gathering more and more varied data. We also hope to expand this to include a wider set of rules covering a larger set of sound changes, and eventually evaluate the method using a human study. We hope that this method of mistake captioning can be applied to other fields of study.


extfieldsofstudy:
Computer Science
Engineering

text:
KUIELab-MDX-Net: A Two-Stream Neural Network for Music Demixing

Recently, many methods based on deep learning have been proposed for music source separation. Some state-of-the-art methods have shown that stacking many layers with many skip connections improve the SDR performance. Although such a deep and complex architecture shows outstanding performance, it usually requires numerous computing resources and time for training and evaluation. This paper proposes a two-stream neural network for music demixing, called KUIELab-MDX-Net, which shows a good balance of performance and required resources. The proposed model has a time-frequency branch and a time-domain branch, where each branch separates stems, respectively. It blends results from two streams to generate the final estimation. KUIELab-MDX-Net took second place on leaderboard A and third place on leaderboard B in the Music Demixing Challenge at ISMIR 2021. This paper also summarizes experimental results on another benchmark, MUSDB18. Our source code is available online.

Introduction
Recently, many methods have been proposed for music source separation. Notably, deep learning approaches Défossez et al. (2021) have become mainstream because of their excellent performance. Some state-of-the-art methods (Choi et al., 2020;Takahashi et al., 2018; have shown that stacking many layers with many skip connections improve the SDR performance. Although a deep and complex architecture shows outstanding performance, it usually requires numerous computing resources and time for training and evaluation. Such disadvantages make them not affordable in a restricted environment where limited resources are provided. For example, some deep models such as LaSAFT-Net  exceed the time limit of the Music Demixing Challenge (MDX) at ISMIR 2021  even if they are the current state of the art on the MUSDB18 (Rafii et al., 2017b) benchmark. This paper presents a source separation model named KUIELab-MDX-Net. We empirically found a good balance of performance and required resources to design KUIElab-MDX-Net. For example, we replaced channel-wise concatenation operations with simple element-wise multiplications for each skip connection between encoder and decoder (i.e., for each U-connection in U-Net). In our prior experiments, it reduced parameters with negligible performance degradation.
Also, we removed the other skip connections, especially, skip connections used in dense blocks (Choi et al., 2020;Takahashi et al., 2018;. We observed that stacked convolutional networks without dense connections followed by Time-Distributed Fully connected layers (TDF) (Choi et al., 2020) could perform comparably to dense blocks without TDFs. TDF, proposed in (Choi et al., 2020), is a sequence of linear layers. It is applied to a given input in the frequency domain to capture frequency-to-frequency dependencies of the target source. Since a single TDF block has the whole receptive field in terms of frequency, injecting TDF blocks into a conventional U-Net (Ronneberger et al., 2015) improves the SDR performance on singing voice separation even with a shallower structure.
By introducing such tricks, we found a computationally efficient and effective model design. As a result, the proposed architecture has a time-frequency branch and a time-domain branch, where each branch separates stems, respectively. It blends results from two streams to generate the final estimation. KUIELab-MDX-Net took second place on leaderboard A and third place on leaderboard B in the Music Demixing Challenge at ISMIR 2021. This paper also summarizes experimental results on another benchmark, MUSDB18.

Frequency Transformation for Source Separation
Some source separation methods (Choi, 2021;Choi et al., 2020;Yin et al., 2020) have adopted Frequency Transformation (FT) to capture frequency-to-frequency dependencies of the target source. Both designed their FT blocks with fully connected layers, also known as linear layers. For example, (Choi et al., 2020) proposed Time-Distributed Fully connected layers (TDF) to capture frequency patterns observed in spectrograms of a singing voice. A TDF block is a sequence of two linear layers. It is applied to a given input in the frequency domain. The first layer downsamples the features to R F/bn , where we denote the number of frequency bins in a given spectrogram feature by F and the bottleneck factor that controls the degree of downsampling by bn.
TFC-TDF-U-Net v1 (Choi et al., 2020) proposed the original TFC-TDF-U-Net for singing voice separation. We call this architecture TFC-TDF-U-Net v1 for the rest of this paper. It adopted a Time-Frequency Convolutions followed by a TDF (TFC-TDF) block as a fundamental building block. By replacing fully connected 2-D convolutional building blocks, conventionally used in U-Net (Ronneberger et al., 2015) with TFC-TDF blocks, it showed a promising performance on singing voice separation tasks of the MUSDB18 (Rafii et al., 2017b) dataset. Also, injecting TDF blocks can enhance separation quality for the other tasks of MUSDB18, as shown in (Choi, 2021). (Choi, 2021) presented how adding TDF blocks improves separation quality by visualizing trained weight matrixes of single-layered TDF blocks (they additionally trained U-Nets with single-layered TDF blocks for weight visualization). As shown in Figure 1, each matrix is trained to analyze timbre features uniquely observed in its instrument by capturing harmonic patterns (i.e., y = α β x). It is also observable that the TDF blocks still performs well on each scale.

Figure 1: Weight matrixes visualization of single-layered TDF blocks
We summarized TFC-TDF-U-Net v1's performance reported in (Choi, 2021) in the experiment section. Since the original TFC-TDF-U-Net v1 is computationally heavy to be evaluated within the time limit of the MDX challenge, we could not submit this, although its performance was promising on the MUSDB18 benchmark. To make an affordable model for the MDX challenge, we empirically found a good balance of performance and required resources. Figure 2 2 depicts the overall flow at inference time: the four U-Net-based separation models (TFC-TDF-U-Net v2) first estimate each source independently, then the Mixer model takes these estimated sources (+ mixture) and outputs enhanced estimated sources. Also, we extract sources with another network based on a time-domain approach, as shown on the right side of Figure 2. We used pretrained Demucs (Défossez et al., 2021) without fine-tuning. Finally, it takes the weighted average for each estimated source, also known as blending (Uhlich et al., 2017). On top of these architectural changes, we also use a different loss function (time-domain l 1 loss) as well as source-specific data preprocessing. As shown in Figure 3, high frequencies above the target source's expected frequency range were cut off from the mixture spectrogram. This way, we can increase n_fft while using the same input spectrogram size (which we needed to constrain for the separation time limit), and using a larger n_fft usually leads to better SDR. It is also why we did not use a multi-target model (a single model that is trained to estimate all four sources), where we could not use source-specific frequency cutting.

Mixer
Although training one separation model for each source can benefit from source-specific preprocessing and model configurations, these models lack the knowledge that they are separating using the same mixture. We thought an additional network that could exploit this knowledge (which we call the Mixer) could further enhance the independently estimated sources. For example, estimated 'vocals' often have drum snare noises left. The Mixer can learn to remove sounds from 'vocals' that are also present in the estimated 'drums' or vice versa.
We only tried very shallow models (such as a single convolution layer) for the Mixer during the MDX Challenge due to the time limit. We look forward to trying more complex models in the KUIELab-MDX-Net: A Two-Stream Neural Network for Music Demixing. Proceedings of the MDX Workshop, 2021. future since even a single 1 × 1 convolution layer was enough to make some improvement on total SDR (Section "Performance on the MUSDB18 Benchmark").

Experimental Results
This section describes the model configurations, STFT parameters, training procedure, and evaluation results on the MUSDB18 benchmark. For training, we used the MUSDB-HQ dataset with default 86/14 train and validation splits.

Configurations and Training
We present a comparison between configurations of TFC-TDF-U-Net v1 and v2 as follows. This applies to all models regardless of the target source (we did not explore different model configurations for each source). In short, v2 is a more shallow but wider model than v1.
# blocks # convs per block bn # freq bins # STFT frames hop size The number of intermediate channels is increased/decreased after down/upsampling layers with a linear factor of 32. Also, as mentioned in Section "TFC-TDF-U-Net v2," we used different n_fft for each source: (6144,4096,16384,8192) for (vocals, drums, bass, other).
All five models (four separation models + Mixer) were optimized with RMSProp with no momentum. We used random chunking and mixing instruments from different songs for data augmentation (Uhlich et al., 2017). We also used data augmentation based on pitch shift and time stretch (Défossez et al., 2021). The overall training procedure can be summarized into two steps: 1. Train single-target separation models (TFC-TDF-U-Net v2) for each source. 2. Train the Mixer while freezing the pretrained weights of the separation models.

Performance on the MUSDB18 Benchmark
We compare our models with current state-of-the-art models on the MUSDB18 benchmark using the SiSEC2018 version of the SDR metric (BSS Eval v4 framewise multi-channel SDR). We report the median SDR over all 50 songs in the MUSDB18 test set. Only models for Leaderboard A were evaluated since our submissions for Leaderboard B uses the MUSDB18 test set as part of the training data.
We summarize the MUSDB18 benchmark performance of KUIELab-MDX-Net. We compare it to recent state-of-the-art models: TFC-TDF-U-Net v1 (Choi et al., 2020), X-UMX (Sawata et al., 2021), Demucs (Défossez et al., 2021), D3Net , ResUNet-Decouple+ (Kong et al., 2021). We also include our baselines to validate our architectural design. Even though our models were downsized for the MDX Challenge, we can see that it gives superior performance over the state-of-the-art models and achieves the best SDR for every instrument except 'bass.' Also, it is notable that TFC-TDF-U-Net v2 with Mixer (i.e., v2 + Mixer) outperforms the existing methods except for 'vocals' even without blending with Demucs.


Black Hole Cannibalism

We consider a version of the Hayden-Preskill thought experiment in which the message thrown into the black hole is itself a smaller black hole. We then discuss the implications of the existence of a recovery channel for this black hole message at asymptotic infinity, resulting in a sharpening of the black hole information paradox for observers who never need to approach a horizon. We suggest decoherence mechanisms as a way of resolving this sharpened paradox.

Introduction
The firewall paradox [1,2] is a tension between four widely accepted postulates about black hole evaporation: 1. Unitary. There is a unitary governing the dynamics of black hole formation and evaporation.
2. Semi-classical field theory. Physics outside the stretched horizon of a large black hole is well approximated by semi-classical field equations of a low energy EFT with local Lorentz invariance.
3. Black Hole Entropy/Discrete Energy Levels. To distant observers, a black hole behaves like a quantum system with discrete energy levels, where the dimension of the subspace of states is the exponential of the Bekenstein-Hawking entropy.

No
Drama. An observer who falls freely across the horizon experiences nothing unusual.
Postulates (1) and (3) imply maximal entanglement between the late and early Hawking radiation, but postulates (2) and (4) imply that a portion of the late Hawking radiation is highly entangled with modes in the black hole interior. The full set of postulates, therefore, implies a violation of monogamy of entanglement for this portion of the late Hawking modes [1]. Proposed resolutions include (i) modifying quantum mechanics [3,4], (ii) allowing a breakdown of no drama [1], (iii) violating unitarity [5], (iv) identifying earlier Hawking radiation with the black hole interior [6], (v) modifying the interior geometry [7][8][9], (vi) invoking quantum complexity theory [10][11][12], (vii) allowing for remnants [13], (viii) violating locality [14,15], and (ix) implementing an auxiliary system to collect the radiation [16,17]. Many of these suggestions-in particular, options (ii), (iv), (v), (vi), and (vii)appeal to the role of the horizon in effectively hiding a portion of the spacetime. In this essay we present a thought experiment that results in a manifest violation of monogamy of entanglement amongst modes at asymptotic infinity, in which none of the modes are hidden behind a horizon. This implies that hiding modes behind a horizon or appealing to strong gravitational effects are not sufficient resolutions. We argue that accounting for decoherence alleviates the tension.

The Thought Experiment
In broad strokes, our thought experiment consists of throwing a small black hole into a larger one and using the Hayden-Preskill protocol to recover the state of the smaller black hole. We argue that, given the above postulates and a qubit model of black holes, this yields a violation of monogamy of entanglement in states held only by the asymptotic observer.
Consider a black hole of mass M 1 formed from the collapse of matter in a pure state. As per Hawking's calculation [18,19], the black hole will radiate maximally entangled pairs, one of which escapes the horizon and the other of which falls into the black hole interior. We model this black hole as a collection of N 1 qubits at its formation and its evaporation as the release of individual qubits into the environment [20]. The interior qubits undergo scrambling dynamics. We can factorize the total Hilbert space as: where H BH 1 (t) corresponds to the Hilbert space of the remaining black hole and H rad 1 (t) denotes that of the Hawking radiation after evolution for some time t. 1 The Hilbert space factors are time-dependent because qubits are transferred from the black hole to the radiation as the black hole evaporates. According to unitarity, the total state |ψ 1 (t) is pure at all times. Let us denote the state of the black hole and its Hawking radiation emitted up to time t as ρ BH 1 (t) and ρ rad 1 (t) respectively, each defined via a partial trace over the neglected system. Assume that the asymptotic observer has been collecting the radiation from this black hole since its formation. Now suppose there is another, smaller black hole in this spacetime with mass M 2 << M 1 , again modeled as a collection of N 2 qubits at its formation, evaporating via release of qubits. Assume that the asymptotic observer collects the radiation from this black hole until, at some time t c after its Page time, it falls into BH 1 . 2 The observer can identify which modes correspond to which black hole by measuring the temperature. 3 After BH 2 has fallen into BH 1 , the asymptotic observer continues collecting radiation emitted by BH 1 . See figure 2 for a summary of this sequence of events.
Hayden and Preskill demonstrated that it is information-theoretically possible for an observer outside a black hole to reconstruct the state of a quantum system (the"message") thrown into the black hole [20]. The observer first collects the Hawking radiation until they possess a system maximally entangled with the black hole, which occurs at the Page time. The message is then thrown into the black hole, which evolves via a deterministic unitary transformation that thoroughly and quickly mixes the message into the black hole's state. As the black hole continues to evaporate, the observer collects the radiation. The outside observer needs to collect only a few more qubits than the length of the original message to reconstruct it.
We choose that the state of the radiation released by BH 1 is maximally entangled with BH 1 at t c . Our asymptotic observer therefore has the ingredients to perform the Hayden-1 Although this evolution is truly time evolution, the choice of time coordinate is arbitrary [12]. In addition, one can use finite dimensional quantum information theory to describe black hole dynamics [21]. 2 One may be concerned that some radiation from BH1 may fall into BH2, preventing it from being collected. We can avoid this by placing the black holes in separate regions connected by an (initially non)traversable wormhole. We can choose the momenta of the black holes and when the wormhole is made traversable such that the black holes collide at the appropriate time. To avoid complications with emission of radiation during transit through the wormhole, we can tune the wormhole to be very short, so that the transit time is much smaller than the expected time elapsed between emission of Hawking quanta. Alternatively, we could construct large Dyson spheres surrounding the black holes, which collect the HR from each of them until their point of collision. 3 Or, the modes can be easily distinguished because the BH's have not yet been sent through the traversable wormholes and are therefore separate. Or if using Dyson spheres, the observer can see which Dyson sphere collected which mode.
Preskill protocol to recover the state of BH 2 when it fell into BH 1 , ρ BH 2 (t c ). Recall that the asymptotic observer also collected the radiation from BH 2 , which we can partition into "early" and "late" radiation corresponding to emission before and after the Page time of BH 2 respectively. Therefore, the asymptotic observer possesses the following states: 1. ρ earlyRad 2 (t c ), containing N 2 /2 + δ qubits and defined from ρ rad 2 (t c ) by tracing out the qubits associated with the late radiation.
2. ρ lateRad 2 (t c ), defined from ρ rad 2 (t c ) by tracing out the qubits associated with the early radiation.
See figure 2 for a circuit-style diagram of this scenario. As we shall see, these systems exhibit a violation of monogamy of entanglement. Page's theorem [22] states that if we bipartition a statistically typical 4 quantum system into A and B with Hilbert space dimensions m and n respectively, and if 1 m ≤ n, then This implies that if the joint system AB is pure and if m and n are very large, then the smaller subsystem is close to maximally mixed. In our setup, the qubits of (2) and (3) are A and the qubits of (1) are B of a pure joint system AB. We can apply Page's theorem by choosing appropriate system sizes. Therefore, the qubits making up subsystems (2) and (3) are maximally entangled with (1). A maximally mixed state has a diagonal density matrix, so ρ 2,3 is roughly diagonal. Tracing out the qubits of (3) thus yields a smaller-dimensional nearly diagonal density matrix for the qubits of (2). Therefore, (2) is nearly maximally entangled with (1). Just before BH 1 swallows BH 2 , Hawking's calculation implies a high degree of entanglement across the horizon of BH 2 . Therefore, the qubits identified with the BH 2 's interior and those identified with its radiation are highly entangled. Crossing the event horizon of BH 1 should not change the entanglement structure of BH 2 's interior. Because the recovery procedure involves only LOCC operations, the entanglement structure of (3) is the same as that of the original state. Therefore, subsystem (3) is highly entangled with (2). Because (2) is maximally entangled with (1), this is a violation of monogamy of entanglement. Therefore, the 3 subsystems held by the asymptotic observer exhibit a violation of monogamy of entanglement.
This thought experiment thus presents a sharp version of the firewall paradox for the asymptotic observer, who is free to manipulate these states, e.g. by encoding them in a spin chain, and has no time constraint on doing so. Therefore, no resolution to this paradox can rely on strong gravitational effects or on hiding something behind the black hole horizon.  Figure 1. The three main steps of the thought experiment. a. The two black holes are initially in separate regions (see footnote 2), emitting Hawking radiation. The asymptotic observer collects the radiation emitted from both black holes. b. The BH 2 has fallen into BH 1 . BH 1 continues to emit radiation, which is collected by the asymptotic observer. The observer runs this radiation (and the radiation previously collected from BH 1 ) through a decoder. c. Finally, all the radiation in the possession of the asymptotic observer (including the reconstructed state of BH 2 ) is encoded into a spin chain.  Figure 2. A diagram of the Hayden-Preskill protocol in which the "message" is another black hole. Black hole 1 (BH 1 ), the larger black hole, has been evaporating for a long time and has become maximally entangled with its previously emitted radiation. The "message" is the state of black hole 2 (BH 2 ), the smaller black hole, which is also maximally entangled with its previously emitted radiation. The message (BH 2 ) falls into BH 1 and the joint system evolves unitarily. More radiation is emitted from BH 1 and collected by the asymptotic observer. The asymptotic observer runs this radiation and the radiation previously emitted from BH 1 through a decoder to recover the state of BH 2 . The asymptotic observer also collects the radiation emitted from BH 2 before it fell into BH 1 .

A Proposed Resolution via Decoherence
In the spirit of [12,25], we argue that accounting for decoherence resolves this paradox. Decoherence is the observation that quantum systems interacting with an environment become entangled with or, equivalently, leak information into said environment [26][27][28][29]. Decoherence causes the system's the density matrix to become nearly diagonal in a special basis, which is determined by details of the interaction. Frequently paired with decoherence is the Everettian view of quantum mechanics, which holds that measurement restricts the observer to a particular "branch" of the wavefunction corresponding to a measurement outcome. 5 The Everettian picture may seem to privilege an arbitrary basis for the decomposition of the global wavefunction into branches, but decoherence provides a particular basis in which the density matrix becomes essentially diagonal, corresponding to a superposition of classical macroscopic states. For evaporating black holes, wavefunction branching occurs when something interacts with its released Hawking quanta, which have, for example, a range of possible outgoing momenta. A "definite geometry" requires a choice of momentum for each emitted Hawking mode and, by extension, a definite position and momentum for the black hole. The total state at non-asymptotic values of t describes an ensemble of these possible geometries [12]. A classical (definite) geometry exists only on a branch of the wavefunction. As the Hawking modes experience interactions, the state of the black hole decoheres into some preferred basis.
Fully defining the process by which BH 2 falls into BH 1 requires a definite point and time at which BH 2 crosses the event horizon of BH 1 , thus requiring a definite geometry, or perhaps a projection to a subset of geometries. Therefore, an instance of this process must occur on a branch (or subset of branches) rather than on the global wavefunction. The asymptotic observer in our setup, however, sees the global wavefunction, which is a superposition of geometries. Therefore, the global observer cannot possess a definite reconstructed state ρ BH 2 (t c ), meaning the global observer does not possess definite states in violation of monogamy of entanglement. Thus, accounting for decoherence resolves the paradox we have outlined in this work. 6 serts that the global wavefunction continues to exist post-measurement, while collapse-based interpretations deny this. Thus in the Everettian view, all measurement possibilities are realized within the global wavefunction as separate branches, hence the alternative name "Many Worlds interpretation." There is actually a family of related views, many of which are called the Everettian interpretation or the Many Worlds interpretation. Sometimes these terms refer to distinct views in this family, and other times they are conflated. The subtleties involved are not critical for our discussion. 6 Our proposed resolution is related to recent works applying the RT formula to a black hole coupled to an auxiliary system [16,17]. The lesson from decoherence for the firewall paradox is that if the bulk observer has access only to the system, their experiencing information loss does not conflict with unitary evolution of the system + environment. The recent works applying the RT formula to a black hole + auxiliary system have the same form: the "environment" is the auxiliary system and the black hole is the "system." Escape of Hawking radiation in [16] with absorbing boundary conditions corresponds to leakage of information into the environment, and the bulk observer has access only to the system. Again, the lesson is that confusing a part (the BH) with the whole (BH + auxiliary system) leads to a contradiction.