extfieldsofstudy:
Medicine
Computer Science

text:
Learning Carbohydrate Digestion and Insulin Absorption Curves Using Blood Glucose Level Prediction and Deep Learning Models

Type 1 diabetes is a chronic disease caused by the inability of the pancreas to produce insulin. Patients suffering type 1 diabetes depend on the appropriate estimation of the units of insulin they have to use in order to keep blood glucose levels in range (considering the calories taken and the physical exercise carried out). In recent years, machine learning models have been developed in order to help type 1 diabetes patients with their blood glucose control. These models tend to receive the insulin units used and the carbohydrate taken as inputs and generate optimal estimations for future blood glucose levels over a prediction horizon. The body glucose kinetics is a complex user-dependent process, and learning patient-specific blood glucose patterns from insulin units and carbohydrate content is a difficult task even for deep learning-based models. This paper proposes a novel mechanism to increase the accuracy of blood glucose predictions from deep learning models based on the estimation of carbohydrate digestion and insulin absorption curves for a particular patient. This manuscript proposes a method to estimate absorption curves by using a simplified model with two parameters which are fitted to each patient by using a genetic algorithm. Using simulated data, the results show the ability of the proposed model to estimate absorption curves with mean absolute errors below 0.1 for normalized fast insulin curves having a maximum value of 1 unit.

Introduction
Type 1 diabetes is among the most frequent chronic diseases in children and young adults and has an important social and health impact that affects the quality of life of the people who suffer from it. The disease begins when the pancreas does not produce enough insulin or when it is not able to use the insulin produced efficiently. As a result, an increase in the level of glucose in the blood can cause heart and kidney problems or damage some of the nerves in the body.
Despite the latest advances, an exact mechanism for curing the disease is not known. This increases the importance of proposing and validating techniques that facilitate the monitoring and treatment of type 1 diabetics [1] and even the development of systems that allow the identification of patterns that can help to find the best treatment for each patient.
Machine learning techniques and algorithms are proving to be able to build a solid foundation in many areas in general and for blood glucose prediction in particular [2,3]. Data from real patients, by applying data-based models and techniques able to learn and extract patterns from the data, have been used for monitoring patients and predicting possible drops or rises in blood glucose levels, establishing recommendations based on each user's particular and personal data. Currently, both physiological inspired models [4] and machine learning models (in particular deep learning-based models [3]) are being successfully applied to data obtained from blood glucose monitors (Continuous Glucose Monitors or CGM) together with manually provided information such as insulin boluses and meals. Physiological-inspired models are mapped to differential equations controlling the blood glucose level kinetics based on user-dependent parameters, which are sometimes hard to estimate based on the data collected from user glucose sensors (such as CGM) and insulin and carbohydrate intakes. Machine learning models represent black-box generic models with some internal parameters that are tuned to fit the data obtained for the sensors in order to minimize the error when using the model for estimating unknown values. Machine learning models try to find patterns in blood glucose levels based on the times, types, and amounts of insulin used and the carbohydrate content taken. Physiological-inspired models [4] show that the human body follows a process that takes an absorption step (both for carbohydrates and insulin) first and then uses the plasma insulin and glucose levels to guide the different insulin and glucose-dependent processes in the human organs and tissues. Intuitively, applying machine learning models directly to insulin bolus units consumed and carbohydrate content in meals taken will make it more difficult to learn blood glucose patterns than using carbohydrate and insulin absorption curves. However, carbohydrate and insulin absorption curves can not be directly measured from the glucose sensors used by type 1 diabetes patients.
This paper proposes a novel mechanism to estimate patient-dependent meal digestion and insulin absorption curves based on data-driven genetic optimization and uses the estimated curves in order to optimize the performance of a Recurrent Neural Network (RNN) based on Long Short-Term Memory (LSTM) cells in order to estimate future levels of blood glucose.

State of the Art
Type 1 diabetes mellitus (T1DM) patients adjust insulin doses in order to properly fit the carbohydrate intake to keep blood glucose levels in range. An optimal regulation of blood glucose levels will require a proper combination of timely injections of insulin boli (considering the number of insulin units and the type of insulin to be used) and food intake depending on the current and past blood glucose levels.
Several absorption models for insulin and carbohydrates have been proposed in order to assess the influence over time of human actions (meal intake or insulin injections) into plasma concentration levels and therefore provide an estimation of the impact on blood glucose levels. Chiara et al. [5] proposed a simulation model in normal humans that describes the physiological events that occur after a meal. The parameters for the model were set to fit the mean data of a large normal subject database that underwent a triple tracer meal protocol that provided quasi-model-independent estimates of major glucose and insulin fluxes. A parametric model was developed decomposing the system into subsystems. Wilinska et al [6] proposed a model for insulin Lispro kinetics with bolus and Continuous Subcutaneous Insulin IInfusion (CSII) modes of insulin delivery. Eleven alternative models of insulin kinetics were proposed. The models considered several absorption approximations, taking into account different rates of absorption and considering the influence of the volumes of insulin. The models also studied two pathways for insulin absorption (fast and slow) and were based on compartment models.
Absorption models for insulin and carbohydrates can be expressed in mathematical terms in the form of differential equations. A particular mathematical representation can be found in Ruan et al [4]. The authors developed a new hierarchical model to relate subcutaneous insulin delivery and carbohydrate intake to continuous glucose monitoring over 12 weeks while describing day-to-day variability. A compartment model that comprised five linear differential equations was proposed. Hovorka et al. [7] also proposed a compartment model, which represented the glucoregulatory system and included submodels representing the absorption of subcutaneously administered short-acting insulin Lispro  [8] presented a Plasma Insulin Concentration (PIC) estimator that incorporated Hovorka's glucose-insulin model [7] with the unscented Kalman filtering algorithm. The authors took an interest in improving convergence. Haiya et al. [9] modeled the insulin therapies using a delay differential equation model. The authors studied the dynamics of the model both qualitatively and quantitatively. Lehmann and Deutsch [10] implemented the equations in the model in [11] for in silico experimentations in the AIDA2 simulator [12]. The authors used numerical integration in order to find solutions to the differential equations guiding the glucose kinetics in a simulated patient. The model contains a single glucose pool representing extracellular glucose (including blood glucose) into which glucose enters via both intestinal absorption and hepatic glucose production [10]. Intestinal absorption is modeled using a first-order equation. Glucose enters the portal circulation using this absorption model from the gut. Insulin absorption is also guided by a differential equation in which the change in plasma insulin concentration is controlled by first-order rate constant elimination and an insulin absorption model that considers the delay in action [10].
Machine learning techniques have also been used as an alternative family of methods in order to estimate upcoming values of blood glucose (BG) levels. Machine learning methods are able to learn from the experience (using labeled training data) how to combine the information measured from glucose, carbohydrate, and insulin signals in order to minimize the errors in the prediction of upcoming blood glucose (BG) levels. Machine learning models are based on multiple parameters that can be trained in order to approximate or simulate the underlying physiological model. However, generic machine learning models are not physiologically related to the human body glucose kinetic processes.
Different machine learning models have been proposed in order to predict future Blood glucose (BG) levels are based on the knowledge of factors such as insulin dosage, nutritional intake, physical exercise, and recent BG levels. Pappada et al. [13] proposed the use of a Neural Network-based model to determine the impact of insulin, meals, daily activities, lifestyle, and emotional states on predicting glycemic trends. The study found that the Neural Network model worked well at predicting normal and hyperglycemic events but provided worse results for hypoglycemia episodes. A different machine learning model based on the implementation of a Support Vector Regression (SVR) algorithm that used a Kalman filter in order to estimate hidden values in the model was proposed in [14]. The model was trained for blood glucose estimation based on current and past levels of carbohydrates, insulin, and BG levels. The authors found that the prediction results provided by the machine learning model were comparable to those manually predicted by a doctor. Another study that uses a machine learning model based on SVR to estimate upcoming BG levels can be found in [15]. The authors limited the input signal to the information provided by a Continuous Glucose Monitoring (CGM) device. The study added a Differential Evolution (DE) algorithm to improve results. A different model based on Artificial Neural Networks (ANNs) for predicting blood glucose levels can be found in [16]. The study also focused on the use of CGM data alone. A different approach for BG estimations using a meta-learning approach can be found in [17]. The results are based on the use of regularized learning algorithms. The paper designed a mechanism that facilitated the portability of the model from patient to patient.
Recently, some deep learning algorithms have also been applied to the prediction of BG levels trying to achieve a better performance compared to previously used methods. The authors in [18] showed some methods for deep multi-output blood glucose forecasting and validated that the results using deep learning methods outperformed previous shallow learning alternatives. Mhaskar et al. [19] also proposed a deep learning approach to BG level estimation based on the previous BG levels but using a pre-clustering mechanism to train specific models for hypo, eu and hyperglycemic segments. The authors also demonstrated that deep learning methods can outperform shallow networks.
In this manuscript, we propose and validate a new mechanism that combines some information from the glucose kinetic physiological processes with the ability of deep machine learning models to adapt to each individual data in order to achieve better glucose level predictions. The model incorporates physiological absorption models for insulin and carbohydrates into the pre-processing steps of a machine learning algorithm using a Recurrent Neural Network based on LSTM cells. Then, the learning algorithm is able to learn the intricate dependencies between blood sugar levels and estimated values for blood insulin and carbohydrate concentrations trying to make it simpler for the deep learning model to learn.

Materials and Methods
This section describes the dataset that has been used in order to evaluate the proposed model. The section also contains the details of the model and the methods used for validation.

Dataset Used
The AIDA2 simulator [10] is based on an in silico model for the glucose-insulin interaction in the human body. The simulator comes with 40 different case scenarios, which can be extended. The simulator generates blood glucose curves as well as plasma insulin curves based on the amount of carbohydrates and insulin doses taken. The AIDA2 simulator uses numerical integration in order to solve the differential equations guiding the glucose kinetics inside a human body [10]. The model is based on 4 differential equations along with 12 auxiliary relations, which are captured in [10]. The derivative of the insulin concentration over time has a positive term proportional to the rate of insulin absorption divided by the volume of insulin and a negative term controlled by the rate of insulin elimination [10]. The change of glucose concentration with time is also guided by a differential equation that takes into account the systemic appearance of glucose via glucose absorption from the gut, the overall rate of peripheral and insulin-independent glucose utilization, the net hepatic glucose balance, the rate of renal glucose excretion, and the volume of distribution for glucose [10].
In this manuscript, we want to analyze if meal digestion and insulin absorption curves can be learned by a proposed machine learning model using genetic optimization. The AIDA2 simulator constitutes a convenient dataset, since it makes the information about insulin absorption curves computed by the simulating model available, which can be used in order to validate the approach proposed in this manuscript. The AIDA2 simulator [12] has also been used in previous research studies such as [20,21] in order to validate blood glucose prediction algorithms. Figure 1 shows the result for a particular simulated patient for the plasma insulin levels from a single simulation in which the levels of plasma insulin are derived from the amounts and types of insulin doses taken (shown in the lower part of the figure). The change in plasma insulin depends on the rate of insulin absorption, the volume of insulin distribution, and a first-order rate constant for insulin elimination. The numerical integration of the model shows plasma insulin curves that have a growing first part in which the insulin absorption is higher than the insulin elimination, and a second decreasing part is controlled by the dominant elimination part.  The characteristics of the dataset used are captured in Table 1. For each case scenario provided by the AIDA 2 simulator [10], a total of 8 days of data is generated. The average number of insulin injections per day has been set to 4, and the average number of meals per day has been set to 6.

Proposed Model
Deep learning models for blood glucose (BG) prediction previously found in the literature such as [18] use the raw information from insulin and carbohydrate intakes directly as inputs to a multilayer model, which tries to estimate the BG levels in a prediction horizon (PH) (normally of 30 or 60 min). The model we propose uses a two-step approach in which insulin and carbohydrate absorption curves are estimated first and then used to feed a multilayer deep learning model.
Insulin and carbohydrate absorption curves can be estimated by solving mathematical models characterizing the glucose metabolism such as [11]. The software in [12] has used numerical integration to solve the differential equations in [11]. Solving the models for each user based on available measured data from Continuous Glucose Monitors (CGM) requires the estimation of internal parameters for each mathematical model, which could be difficult in some cases. In order to simplify the machine learning process to fit the available data, we propose an absorption curve generation process that will use a linear approximation for the first ascending part of the curve as used in [9] and an exponential decay for the second descending part of the curve as in [22].
A deep learning Recurrent Neural Network (RNN) model based on LSTM cells [19] will be used in order to generate estimates for the upcoming levels for blood glucose (BG) based on current and past levels for the same BG signal and the output of the 3 curve generation processes, one for the digestion of carbohydrates and the other two for fast and slow insulin absorption. The curve generators will receive as parameters the time to peak and the after-peak decay rate in order to generate signals, which will approximate the The characteristics of the dataset used are captured in Table 1. For each case scenario provided by the AIDA 2 simulator [10], a total of 8 days of data is generated. The average number of insulin injections per day has been set to 4, and the average number of meals per day has been set to 6.

Proposed Model
Deep learning models for blood glucose (BG) prediction previously found in the literature such as [18] use the raw information from insulin and carbohydrate intakes directly as inputs to a multilayer model, which tries to estimate the BG levels in a prediction horizon (PH) (normally of 30 or 60 min). The model we propose uses a two-step approach in which insulin and carbohydrate absorption curves are estimated first and then used to feed a multilayer deep learning model.
Insulin and carbohydrate absorption curves can be estimated by solving mathematical models characterizing the glucose metabolism such as [11]. The software in [12] has used numerical integration to solve the differential equations in [11]. Solving the models for each user based on available measured data from Continuous Glucose Monitors (CGM) requires the estimation of internal parameters for each mathematical model, which could be difficult in some cases. In order to simplify the machine learning process to fit the available data, we propose an absorption curve generation process that will use a linear approximation for the first ascending part of the curve as used in [9] and an exponential decay for the second descending part of the curve as in [22].
A deep learning Recurrent Neural Network (RNN) model based on LSTM cells [19] will be used in order to generate estimates for the upcoming levels for blood glucose (BG) based on current and past levels for the same BG signal and the output of the 3 curve generation processes, one for the digestion of carbohydrates and the other two for fast and slow insulin absorption. The curve generators will receive as parameters the time to peak and the after-peak decay rate in order to generate signals, which will approximate  Table 2. human body absorption processes. The model is shown in Figure 2. The details for the RNN part of the model are captured in Table 2. The AIDA2 simulator provides data samples every 15 min. In order to feed the model in Figure 2, the input data are pre-processed in order to generate 4-h windows of data for each of the input signals. No further pre-processing steps have been required, since the data contain no measurement errors. The output of the model in Figure 2 will be the estimated level for blood glucose in a prediction horizon (PH) of 60 min. The RNN part of the model in Figure 2 will receive the generated absorption curves together with the values measured from the CGM device as the input information in order to learn patterns in the data and fit the weights to minimize the output estimation error. The RNN model is expected to be able to reduce the estimated errors for blood glucose levels in the PH when the generated input signals optimally fit the real carbohydrate and insulin blood concentrations, since these values are the ones used in the internal model used to generate the real output data. In order to optimally learn the absorption curves for carbohydrates and slow and fast insulin, the RNN model in Figure 2 will be trained using the output for each of the 3 curve generator processes using different values of peak times and decay rates. In order to generate the 3 absorption curves, a total of 6 parameters (3 peak times and 3 decay rates) have to be set. The RNN part of the model in Figure 2 will then be trained, and the output error after training will be used as an accuracy value for validating the similarity between the generated curves and the real ones. A genetic algorithm is used in order to find the optimal values for the 6 parameters controlling the generation of the absorption curves. In The AIDA2 simulator provides data samples every 15 min. In order to feed the model in Figure 2, the input data are pre-processed in order to generate 4-h windows of data for each of the input signals. No further pre-processing steps have been required, since the data contain no measurement errors. The output of the model in Figure 2 will be the estimated level for blood glucose in a prediction horizon (PH) of 60 min. The RNN part of the model in Figure 2 will receive the generated absorption curves together with the values measured from the CGM device as the input information in order to learn patterns in the data and fit the weights to minimize the output estimation error. The RNN model is expected to be able to reduce the estimated errors for blood glucose levels in the PH when the generated input signals optimally fit the real carbohydrate and insulin blood concentrations, since these values are the ones used in the internal model used to generate the real output data. In order to optimally learn the absorption curves for carbohydrates and slow and fast insulin, the RNN model in Figure 2 will be trained using the output for each of the 3 curve generator processes using different values of peak times and decay rates. In order to generate the 3 absorption curves, a total of 6 parameters (3 peak times and 3 decay rates) have to be set. The RNN part of the model in Figure 2 will then be trained, and the output error after training will be used as an accuracy value for validating the similarity between the generated curves and the real ones. A genetic algorithm is used in order to find the optimal values for the 6 parameters controlling the generation of the absorption curves. In particular, the algorithm described in chapter 7 in reference [23] is used. The algorithm concatenates 4 major functions in order to find the individuals that optimize the value of a fitness function: recombine, mutate, evaluate, and select. The recombine function generates individuals that combine the attributes from those in the previous generation in order to find new individuals. The mutate function makes some random changes with a particular probability to the attributes in each individual. The evaluate function calculates the fitness values for all the individuals in a population. Finally, the selection function chooses the best individuals in order to generate a new population. The algorithm can iterate several times or epochs. An initial population of 15 individuals has been used, including a mutation probability of 20% in order to generate 15 generations of individuals. The best-fitted individuals have been selected for each generation. The process is shown in Figure 3. The initial values for each individual have been randomly selected. The maximum and minimum values for peak times have been set from 0 to 2 h for meal digestion and fast insulin concentration and from 0 to 12 h for slow insulin absorption.
Sensors 2021, 21, x FOR PEER REVIEW 7 of 14 particular, the algorithm described in chapter 7 in reference [23] is used. The algorithm concatenates 4 major functions in order to find the individuals that optimize the value of a fitness function: recombine, mutate, evaluate, and select. The recombine function generates individuals that combine the attributes from those in the previous generation in order to find new individuals. The mutate function makes some random changes with a particular probability to the attributes in each individual. The evaluate function calculates the fitness values for all the individuals in a population. Finally, the selection function chooses the best individuals in order to generate a new population. The algorithm can iterate several times or epochs. An initial population of 15 individuals has been used, including a mutation probability of 20% in order to generate 15 generations of individuals. The bestfitted individuals have been selected for each generation. The process is shown in Figure  3. The initial values for each individual have been randomly selected. The maximum and minimum values for peak times have been set from 0 to 2 h for meal digestion and fast insulin concentration and from 0 to 12 h for slow insulin absorption.

Software Configuration
The model in Figures 2 and 3 has been implemented in Python. The DEAP library [24] has been used in order to implement the genetic algorithm-based optimization process for estimating the optimal parameters used to generate the absorption curves. Each individual in the population is defined by the values for the six parameters controlling the generation of the curves for insulin and carbohydrate absorption processes (peak times and decay rates). A FitnessMax function has been used with weights equals to -1 so that a minimization problem is solved. An initial population of 15 individuals has been selected. The model in Figure 2 is used as the fitness function so that the individuals able to generate curves which are able to better train the RNN model are preferred. The cxTwoPoint parameter is used to cross individuals. The mutGaussian parameter is used to define the mutation of the individuals. A mutation probability of 20% has been selected.
The parameters for the LSTM-based Recurrent Neural Network in Figure 2 are captured in Figure 4. Two LSTM layers are stacked in order to augment the capacity of the model to learn intricate patterns and dependencies from the input data. Two-layer models have been previously used in studies such as [19], showing good results. The number of memory units in the LSTM cells for the first layer has been set to 10 and has been reduced to 5 in the second LSTM layer. These values have been previously validated in studies

Software Configuration
The model in Figures 2 and 3 has been implemented in Python. The DEAP library [24] has been used in order to implement the genetic algorithm-based optimization process for estimating the optimal parameters used to generate the absorption curves. Each individual in the population is defined by the values for the six parameters controlling the generation of the curves for insulin and carbohydrate absorption processes (peak times and decay rates). A FitnessMax function has been used with weights equals to -1 so that a minimization problem is solved. An initial population of 15 individuals has been selected. The model in Figure 2 is used as the fitness function so that the individuals able to generate curves which are able to better train the RNN model are preferred. The cxTwoPoint parameter is used to cross individuals. The mutGaussian parameter is used to define the mutation of the individuals. A mutation probability of 20% has been selected.
The parameters for the LSTM-based Recurrent Neural Network in Figure 2 are captured in Figure 4. Two LSTM layers are stacked in order to augment the capacity of the model to learn intricate patterns and dependencies from the input data. Two-layer models have been previously used in studies such as [19], showing good results. The number of memory units in the LSTM cells for the first layer has been set to 10 and has been reduced to 5 in the second LSTM layer. These values have been previously validated in studies such as [21]. Dropout regularization layers have also been used in order to minimize overfitting problems. Two dense layers are connected at the output of the second LSTM layer in order to generate a single prediction value to estimate the level of blood glucose in the prediction horizon of 60 min. A Keras sequential model [25] has been used in order to implement the model in Figure 4. such as [21]. Dropout regularization layers have also been used in order to minimize overfitting problems. Two dense layers are connected at the output of the second LSTM layer in order to generate a single prediction value to estimate the level of blood glucose in the prediction horizon of 60 min. A Keras sequential model [25] has been used in order to implement the model in Figure 4.

Numerical Results
The AIDA2 simulator [10] has been used to simulate data for 8 consecutive days for 40 different reference T1DM patients. The model in Figures 2 and 3 has been used in order to estimate the absorption curves for each particular patient. The parameters defining the differential equations implemented by the AIDA2 simulator [10] for each patient consider different characteristics that try to accommodate real cases for different users. Therefore, the model in Figures 2 and 3 has been used for each patient independently. The AIDA2 simulator [10] provides access not only to the time evolution of blood glucose levels but also the plasma insulin concentration levels over time. We have used the insulin curves provided by the simulator as the desired or optimal curves to feed the RNN network in our model, since these curves are the real ones used by the simulator in order to calculate future blood glucose (BG) levels. The mean absolute error between the plasma insulin curve provided by the simulator and the approximated curve learned by the proposed model has been computed for each simulated use case. The average results for the learned insulin curves are captured in Table 3. The graphical representation for the AIDA2 provided and the learned curves are shown in Figure 5. Both curves have been normalized so that their maximum value is 1. The results in Figure 5 show that the genetic algorithm obtains a generated insulin curve that has a peak 90 min after the bolus injection and the

Numerical Results
The AIDA2 simulator [10] has been used to simulate data for 8 consecutive days for 40 different reference T1DM patients. The model in Figures 2 and 3 has been used in order to estimate the absorption curves for each particular patient. The parameters defining the differential equations implemented by the AIDA2 simulator [10] for each patient consider different characteristics that try to accommodate real cases for different users. Therefore, the model in Figures 2 and 3 has been used for each patient independently. The AIDA2 simulator [10] provides access not only to the time evolution of blood glucose levels but also the plasma insulin concentration levels over time. We have used the insulin curves provided by the simulator as the desired or optimal curves to feed the RNN network in our model, since these curves are the real ones used by the simulator in order to calculate future blood glucose (BG) levels. The mean absolute error between the plasma insulin curve provided by the simulator and the approximated curve learned by the proposed model has been computed for each simulated use case. The average results for the learned insulin curves are captured in Table 3. The graphical representation for the AIDA2 provided and the learned curves are shown in Figure 5. Both curves have been normalized so that their maximum value is 1. The results in Figure 5 show that the genetic algorithm obtains a generated insulin curve that has a peak 90 min after the bolus injection and the decay rate is also learned to mimic the behavior of the curve used by the simulator (Figure 5 shows that both curves coincide around 3 and a half h after the bolus injection). decay rate is also learned to mimic the behavior of the curve used by the simulator ( Figure  5 shows that both curves coincide around 3 and a half h after the bolus injection).  The genetic algorithm converges after 10 iterations (as shown in Figure 6). Each individual in the population contains the peak values and decay rates, which are used to generate the insulin and carbohydrate absorption curves, which together with the past blood glucose (BG) values are used as the input data to feed the RNN model in Figure 2. The input data are divided into training and validation subsets (70%/30%). The training samples are used to train the RNN model in Figure 2. The validation subset is used to compute the output value for the fitness function that the genetic algorithm uses in the optimization process. The initial population (iteration 0 in Figure 6) is composed of randomly generated values for the input parameters in order to generate random initial absorption curves. After training the RNN in Figure 2 with the different generated curves in the initial population, a significant difference is observed among the different prediction errors (mean value of 1, min value of 0.55, and max value of 2.45). Therefore, the shapes of input curves are important in order for the RNN to be able to predict upcoming values of blood glucose levels based on the past values of insulin and carbohydrate intakes. The best individuals are crossed, and a 20% mutation is introduced in order to generate a new generation of individuals for the next epoch in the optimization process performed by the genetic algorithm. Figure 6 captures the output values for the population used in each of the epochs for one particular patient generated using the AIDA2 simulator. After 10 iterations, the genetic algorithm is not able to significantly improve the accuracy achieved by the RNN in the model. The genetic algorithm converges after 10 iterations (as shown in Figure 6). Each individual in the population contains the peak values and decay rates, which are used to generate the insulin and carbohydrate absorption curves, which together with the past blood glucose (BG) values are used as the input data to feed the RNN model in Figure 2. The input data are divided into training and validation subsets (70%/30%). The training samples are used to train the RNN model in Figure 2. The validation subset is used to compute the output value for the fitness function that the genetic algorithm uses in the optimization process. The initial population (iteration 0 in Figure 6) is composed of randomly generated values for the input parameters in order to generate random initial absorption curves. After training the RNN in Figure 2 with the different generated curves in the initial population, a significant difference is observed among the different prediction errors (mean value of 1, min value of 0.55, and max value of 2.45). Therefore, the shapes of input curves are important in order for the RNN to be able to predict upcoming values of blood glucose levels based on the past values of insulin and carbohydrate intakes. The best individuals are crossed, and a 20% mutation is introduced in order to generate a new generation of individuals for the next epoch in the optimization process performed by the genetic algorithm. Figure 6 captures the output values for the population used in each of the epochs for one particular patient generated using the AIDA2 simulator. After 10 iterations, the genetic algorithm is not able to significantly improve the accuracy achieved by the RNN in the model. The patient in Figures 5 and 6 uses a combination of slow and fast-acting insulin. In order to assess the fast insulin curve learning process minimizing the effect of other types of insulin in the data, a patient using only fast-acting insulin has been simulated. The results are captured in Figure 7. Bolus injection takes place at instant 0 in the curve. The linear approximation for the first part of the curve follows the real curve, achieving a maximum value at the same time. The decay rate in this case is a bit faster than the real curve used by the simulator. Complex machine learning models using several LSTM layers have shown better performance when predicting upcoming values of blood glucose (BG) levels than shallow models [3,18,19,21]. A good machine learning model will be able to compensate some inaccuracy levels in the input data. In order to assess the influence of the use of simpler machine learning models for the performance of the genetic algorithm in finding absorption curves that try to fit the real curves used by the AIDA2 simulator, an RNN model The patient in Figures 5 and 6 uses a combination of slow and fast-acting insulin. In order to assess the fast insulin curve learning process minimizing the effect of other types of insulin in the data, a patient using only fast-acting insulin has been simulated. The results are captured in Figure 7. Bolus injection takes place at instant 0 in the curve. The linear approximation for the first part of the curve follows the real curve, achieving a maximum value at the same time. The decay rate in this case is a bit faster than the real curve used by the simulator. The patient in Figures 5 and 6 uses a combination of slow and fast-acting insulin. In order to assess the fast insulin curve learning process minimizing the effect of other types of insulin in the data, a patient using only fast-acting insulin has been simulated. The results are captured in Figure 7. Bolus injection takes place at instant 0 in the curve. The linear approximation for the first part of the curve follows the real curve, achieving a maximum value at the same time. The decay rate in this case is a bit faster than the real curve used by the simulator. Complex machine learning models using several LSTM layers have shown better performance when predicting upcoming values of blood glucose (BG) levels than shallow models [3,18,19,21]. A good machine learning model will be able to compensate some inaccuracy levels in the input data. In order to assess the influence of the use of simpler machine learning models for the performance of the genetic algorithm in finding absorption curves that try to fit the real curves used by the AIDA2 simulator, an RNN model Complex machine learning models using several LSTM layers have shown better performance when predicting upcoming values of blood glucose (BG) levels than shallow models [3,18,19,21]. A good machine learning model will be able to compensate some inaccuracy levels in the input data. In order to assess the influence of the use of simpler machine learning models for the performance of the genetic algorithm in finding absorption curves that try to fit the real curves used by the AIDA2 simulator, an RNN model with a single LSTM layer with only five memory cells has been used. The results are captured in the next subsection.

Using a Simple Model Results
The model in Figure 8 consisting of an RNN with a single LSTM layer has been used instead of the model in Figure 4 in order to assess the influence of the machine learning model when trying to find optimal absorption curves. Previous studies such as [3,18,19,21] have shown that deeper machine learning models are able to increase the accuracy of blood glucose level predictions, both for horizons of 30 and 60 min, compared to models with a single LSTM layer such as the one in Figure 8. However, the use of a simpler machine learning model with a fewer number of parameters to tune will require less memory and training time and will therefore reduce the time needed to find the optimal absorption curves. The model in Figure 8 uses a single LSTM layer with only five memory cells and a final single neuron output layer to generate the blood glucose level prediction. A dropout layer is added as a regularization layer to the model. with a single LSTM layer with only five memory cells has been used. The results are captured in the next subsection.

Using a Simple Model Results
The model in Figure 8 consisting of an RNN with a single LSTM layer has been used instead of the model in Figure 4 in order to assess the influence of the machine learning model when trying to find optimal absorption curves. Previous studies such as [3,18,19,21] have shown that deeper machine learning models are able to increase the accuracy of blood glucose level predictions, both for horizons of 30 and 60 min, compared to models with a single LSTM layer such as the one in Figure 8. However, the use of a simpler machine learning model with a fewer number of parameters to tune will require less memory and training time and will therefore reduce the time needed to find the optimal absorption curves. The model in Figure 8 uses a single LSTM layer with only five memory cells and a final single neuron output layer to generate the blood glucose level prediction. A dropout layer is added as a regularization layer to the model. The results in Figure 9 show the validation error for the first 10 epochs for the same patient generated by the AIDA2 simulator previously shown in Figure 5. The genetic algorithm is able to converge in 10 epochs using the simpler model. However, the average validation error in this case is 0.29 mmol/L (5.22 mg/dL), which is significantly higher than the validation error of 0.095 mmol/L (1.71 mg/dL) achieved by the two-layer model presented in Section 4.1. This result confirms previous results in studies such as [3,18,19,21] in which the impact of the model for predicting upcoming values of blood glucose levels is assessed. Comparing the prediction results with those in Table 4 in reference [21], the addition of absorption curves fitted to each patient, as proposed in this paper, is able to improve the Root Mean Square Error (RMSE) when using similar LSTM-based deep learning models. The authors in [21] achieved an RMSE value of 3.45 mg/dL using the same dataset as in this paper but only for a prediction horizon of 30 min and an optimized deep learning model. Using the two-layer model in Section 4.1, the absorption curves are able to improve the RMSE value to 1.71 mg/dL for a prediction horizon of 60 min.
In order to estimate the quality of the absorption curves found by the genetic algorithm when using the simpler RNN model, the same comparison with the real fast insulin absorption curves used by the AIDA2 simulator has been performed. The results are captured in Table 4. The average time to peak found is now 1.65 h, which is a bit higher than the 1.52 h found when using the two-layer RNN. Moving the peak time apart from the bolus injection is consistent with the estimated average value for the decay rate, which is The results in Figure 9 show the validation error for the first 10 epochs for the same patient generated by the AIDA2 simulator previously shown in Figure 5. The genetic algorithm is able to converge in 10 epochs using the simpler model. However, the average validation error in this case is 0.29 mmol/L (5.22 mg/dL), which is significantly higher than the validation error of 0.095 mmol/L (1.71 mg/dL) achieved by the two-layer model presented in Section 4.1. This result confirms previous results in studies such as [3,18,19,21] in which the impact of the model for predicting upcoming values of blood glucose levels is assessed. Comparing the prediction results with those in Table 4 in reference [21], the addition of absorption curves fitted to each patient, as proposed in this paper, is able to improve the Root Mean Square Error (RMSE) when using similar LSTM-based deep learning models. The authors in [21] achieved an RMSE value of 3.45 mg/dL using the same dataset as in this paper but only for a prediction horizon of 30 min and an optimized deep learning model. Using the two-layer model in Section 4.1, the absorption curves are able to improve the RMSE value to 1.71 mg/dL for a prediction horizon of 60 min.
In order to estimate the quality of the absorption curves found by the genetic algorithm when using the simpler RNN model, the same comparison with the real fast insulin absorption curves used by the AIDA2 simulator has been performed. The results are captured in Table 4. The average time to peak found is now 1.65 h, which is a bit higher than the 1.52 h found when using the two-layer RNN. Moving the peak time apart from the bolus injection is consistent with the estimated average value for the decay rate, which is also higher than the decay rate found for the two-layer RNN (0.77 h −1 vs. 0.57 h −1 ). The quality of the generated absorption curve for the fast insulin is measured by computing the mean absolute error for the difference between the generated curve and the curve computed by the AIDA2 simulator. Table 4 captures a value of 0.088 for the single layer RNN model which is higher than the 0.077 value in Table 3 for the two-layer RNN model. also higher than the decay rate found for the two-layer RNN (0.77 h −1 vs. 0.57 h −1 ). The quality of the generated absorption curve for the fast insulin is measured by computing the mean absolute error for the difference between the generated curve and the curve computed by the AIDA2 simulator. Table 4 captures a value of 0.088 for the single layer RNN model which is higher than the 0.077 value in Table 3 for the two-layer RNN model.   Figure 10 shows the real vs. the estimated absorption curves when using the singlelayer RNN model using the same patient as in Figure 5. The use of the simpler model is able to find the same peak time, but the values estimated for the decay rate are higher than the value in Figure 5, and the mean absolute error therefore increases.    Figure 10 shows the real vs. the estimated absorption curves when using the singlelayer RNN model using the same patient as in Figure 5. The use of the simpler model is able to find the same peak time, but the values estimated for the decay rate are higher than the value in Figure 5, and the mean absolute error therefore increases. also higher than the decay rate found for the two-layer RNN (0.77 h −1 vs. 0.57 h −1 ). The quality of the generated absorption curve for the fast insulin is measured by computing the mean absolute error for the difference between the generated curve and the curve computed by the AIDA2 simulator. Table 4 captures a value of 0.088 for the single layer RNN model which is higher than the 0.077 value in Table 3 for the two-layer RNN model.   Figure 10 shows the real vs. the estimated absorption curves when using the singlelayer RNN model using the same patient as in Figure 5. The use of the simpler model is able to find the same peak time, but the values estimated for the decay rate are higher than the value in Figure 5, and the mean absolute error therefore increases.

Conclusions
This manuscript has presented a method to estimate insulin and carbohydrate absorption curves based on the ability of an RNN machine learning model to achieve optimal Sensors 2021, 21, 4926 13 of 14 predictions for upcoming values of blood glucose levels. Different absorption curves are generated and fed into the RNN in order to select the curves with better blood glucose level predictions over a 60-min horizon. A simple curve generation method using the concatenation of an increasing segment followed by a decreasing segment has been used based on previous studies. The generated curves perform a linear approximation for the increasing blood carbohydrate concentration segment after a meal or the plasma insulin level after a bolus injection. The decaying phase, after the peak value is reached, is approximated using an exponential curve. The peak instant of time and the exponential decay rate are parameters for the curve generation model. The paper has validated the use of a genetic optimization process in order to find the optimal values for the parameters generating the absorption curves. The validation scores for two different RNN models for predicting the blood glucose level in a 60-min prediction horizon have been used as the fitness function for the genetic optimization.
Using a two-layer LSTM-based RNN model provides a better estimation for the fast insulin absorption curves than using a single-layer RNN model. The optimal curves for fast-acting insulin found by the genetic algorithm when using a two-layer RNN are able to estimate the peak instant of time with less than 10 min error. The mean absolute error between the curves provided by the AIDA2 simulator and the optimal curves found by the genetic algorithm is below 0.1 for curves normalized so that the maximum value is 1. Using the two-layer RNN provides an average mean absolute error of 0.077, while the single layer RNN model generates an average mean absolute error of 0.088. This manuscript has used the data produced by the AIDA2 simulator, since the simulator provides information about the plasma insulin levels computed for each scenario. These plasma insulin curves made available by the simulator have been used in order to validate the results obtained by the proposed method in this manuscript. As a future study, a real patient dataset will be used in order to assess the improvement in blood glucose estimations when using a prior computation for carbohydrate and insulin absorption curves as proposed in this manuscript.


Knowledge Transfer by Discriminative Pre-training for Academic Performance Prediction

The needs for precisely estimating a student's academic performance have been emphasized with an increasing amount of attention paid to Intelligent Tutoring System (ITS). However, since labels for academic performance, such as test scores, are collected from outside of ITS, obtaining the labels is costly, leading to label-scarcity problem which brings challenge in taking machine learning approaches for academic performance prediction. To this end, inspired by the recent advancement of pre-training method in natural language processing community, we propose DPA, a transfer learning framework with Discriminative Pre-training tasks for Academic performance prediction. DPA pre-trains two models, a generator and a discriminator, and fine-tunes the discriminator on academic performance prediction. In DPA's pre-training phase, a sequence of interactions where some tokens are masked is provided to the generator which is trained to reconstruct the original sequence. Then, the discriminator takes an interaction sequence where the masked tokens are replaced by the generator's outputs, and is trained to predict the originalities of all tokens in the sequence. Compared to the previous state-of-the-art generative pre-training method, DPA is more sample efficient, leading to fast convergence to lower academic performance prediction error. We conduct extensive experimental studies on a real-world dataset obtained from a multi-platform ITS application and show that DPA outperforms the previous state-of-the-art generative pre-training method with a reduction of 4.05% in mean absolute error and more robust to increased label-scarcity.

INTRODUCTION
Predicting a student's future academic performance is a fundamental task for developing modern Intelligent Tutoring  Figure 1: Interactive features, such as student response and elapsed time for the response, are automatically recorded to the database whenever a student interacts with ITS. On the other hand, more complicated steps are necessary to obtain a test score: a student should take the test in the designated test center, receive the test score, and report the score to ITS. System (ITS) which aims to provide personalized learning experience by supporting educational needs of each individual. However, labels for academic performance, such as test scores, are often scarce since they are external to ITS. For example, as shown in Figure 1, test scores are not automatically collected inside of ITS. Obtaining a test score requires a student to take the test in the designated test center, receive the score, and report the score to ITS. Transfer learning is a commonly taken approach to address such label-scarcity problems across different domains of machine learning. In this framework, a model is first pre-trained to optimize auxiliary objectives with abundant data, and then fine-tuned on the task of interest. In Artificial Intelligence in Education (AIEd) community, [3] introduced Assessment Modeling (AM), a set of pre-training tasks for label-scarce educational problems including academic performance prediction. AM proposed a pre-training method where first, a masked interaction sequence is generated by replacing a set of interactive features which can serve as criteria for pedagogi-cal evaluation with artificial mask tokens. Then, given the masked interaction sequence, a model is pre-trained to predict the masked interactive features. The idea was borrowed from the Masked Language Modeling (MLM) pre-training method proposed in [7]. In the MLM pre-training method, given a masked word sequence where some words in the sequence are replaced with an artificial mask token, a model is pre-trained to predict the masked words. However, recently, [6] pointed out that the MLM pre-training method has poor sample efficiency and suffers from pre-train/fine-tune discrepancy due to the artificial mask token, and proposed a new discriminative pre-training method. Considering the problems are also inherent in AM, potential gains are expected to be obtainable when the discriminative pre-training method is applied to academic performance prediction.

. R e p o r t a s c o r e A u t o m a t i c a l l y c o l l e c t e d
To this end, we propose DPA, a transfer learning framework with Discriminative Pre-training tasks for Academic performance prediction. There are two models in DPA: a generator and a discriminator. In DPA's pre-training phase, the generator is trained to predict the masked interactive features in the same way as AM. Then, given a replaced interaction sequence which is generated by replacing the masked features with the generator's outputs, the discriminator is trained to predict whether each token in the sequence is the same as the one in the original interaction sequence. After the pre-training, the generator is thrown away and only the discriminator is fine-tuned on academic performance prediction. Also, we investigate diverse pre-training tasks for the generator and show that pre-training the generator to predict a student's response is more effective than to predict the correctness and timeliness of their response which were considered as the most pedagogical interactive features in AM. Extensive experimental studies conducted on a realworld dataset collected from a multi-platform ITS application show that DPA outperforms AM with a reduction of 4.05% in mean absolute error and more robust when the degree of label-scarcity increases. Through a series of ablation experiments, we show that DPA's sample efficient pretraining contributed most to the improvement from AM to DPA.

SANTA: A SELF-STUDY SOLUTION EQUIPPED WITH AN AI TUTOR FOR ENGLISH EDUCATION
In this paper, we conduct experiments on a real-world dataset obtained from Santa 1 , a multi-platform ITS with more than a million users in South Korea available through Android, iOS, and Web that exclusively focuses on the Test of English for International Communication (TOEIC) standardized examination. The publicly accessible version of the dataset was released under the name EdNet [4]. The TOEIC consists of two timed sections, Listening Comprehension (LC) and Reading Comprehension (RC). There are a total of 100 multiple choice exercises in each section, and the total score for each section is 495 in steps of 5 points. Santa provides learning experiences of solving exercises, studying explanations, and watching lectures. When a student consumes a specific learning content, Santa diagnoses their current academic status based on their learning activities records and recommends another learning content appropriate for their 1 https://aitutorsanta.com current position. Santa records diverse types of interactive features, such as student response, the duration of time the student took to respond, and the time interval between the current and previous learning activities. However, unlike the interactive features automatically collected from Santa, obtaining the official TOEIC score requires more steps: a student should register and pay for the test, take the test in the designated test center, receive the test score from the Educational Testing Service, and report the score to Santa (Figure 1). Santa collected students' TOEIC score data by offering small gifts to students when they report their scores.

TRANSFER LEARNING FOR ACADEMIC TEST PERFORMANCE PREDICTION
To overcome the label-scarcity problem in academic test performance prediction, we consider burgeoning machine learning discipline of transfer learning. There is an open issue of what information to transfer or which pre-training task is the most effective for academic test performance prediction. Previous studies proposed two types of pre-training methods for AIEd Tasks: interaction-based method which models students' dynamic learning behaviors [13,16,8,3], and content-based method which learns representations of learning contents [15,25,21,26,32]. [3] showed that interactionbased pre-training method outperforms content-based pretraining methods when the pre-trained model is fine-tuned on several label-scarce educational tasks including academic test performance prediction. Following this line of research, we propose a transfer learning framework where a model is pre-trained using only student interaction data, and finetune the pre-trained model on academic test performance prediction. In this paper, we consider the following interactive features: • eid : A unique ID assigned to an exercise solved by a student. There are a total of 14419 exercises in the dataset.
• part: Each exercise belongs to a specific part that represents the type of the exercise. There are a total of 7 parts in the TOEIC.
• response: Since the TOEIC consists of multiple choice exercises and there are four options for each exercise, a student response for a given exercise is one of the options, 'a', 'b', 'c', or 'd'.
• correctness: Whether a student responded correctly to a given exercise. Note that correctness is a coarse version of response since correctness is processed by comparing response with a correct answer for a given exercise.
• elapsed time: The amount of time a student spent on solving a given exercise.
• timeliness: Whether a student responded to a given exercise under the time limit. Note that timeliness is a coarse version of elapsed time since timeliness is processed by comparing elapsed time with the time limit recommended by domain experts for a given exercise.
• exp time: The amount of time a student spent on studying an explanation for an exercise they had solved.  Figure 2: The overall pre-training/fine-tuning process of DPA when each token in an interaction sequence is a set of eid, part, and response, and response is a feature being masked. mask and cls are special tokens for mask and classification, respectively, which are the same as the ones used in [7].
• inactive time: The time interval between the current and previous interactions.
In our experiments, we normalize the values of elapsed time, exp time, and inactive time so they are between 0 and 1 to stabilize the training process.  Figure 2, since the generator's outputs for the masked features are 'b' and 'a', a replaced interaction sequence is [(e419, part4, b), (e23, part3, b), (e4324, part3, a), (e5233, part1, a)]. Then, the discriminator takes the replaced interaction sequence I R as an input, and predicts whether each token in the sequence is the same as the one in the original interaction sequence (original) or not (replaced). After the pre-training, we throw away the generator and fine-tune the pre-trained discriminator on academic test performance prediction. We provide detailed explanations of each component in the generator and the discriminator, and training objective functions in the following subsections.

Interaction Embeddings
The embedding layer produces a sequence of interaction embedding vectors by mapping each interactive feature to an appropriate embedding vector. We take two different approaches to embed the interactive features depending on whether they are categorical (eid, part, response, correctness, and timeliness) or continuous (elapsed time, exp time, and inactive time) variables. If an interactive feature is a categorical variable, we assign unique latent vectors to possible values of the feature including special values for mask (mask ) and classification (cls). Take response as an example, there is an embedding matrix Eresponse ∈ R 6×d emb where each row vector is assigned to one of 'a', 'b', 'c', 'd', mask, and cls. If an interactive feature is a continuous variable, we assign a single latent vector to the feature. Then, an embedding vector for the feature is computed by multiplying the latent vector and a value of the feature. For instance, we compute an embedding vector for elapsed time as et * E elapsed time , where et is a specific value and E elapsed time ∈ R d emb is a latent vector assigned to elapsed time. Also, mask and classification for the continuous interactive features are indicated by setting their values to -1 and 0, respectively. Not only embeddings for interactive features, positional embeddings are also incorporated into Transformer-based models [29] to consider chronological order of each token. Rather than using conventional positional embeddings which stores an embedding vector for every possible position, we adopt axial positional embeddings [19] to further reduce memory usage. The final interaction embedding vector of dimension d emb for each time-step is the sum of all embedding vectors in the time-step. The interaction embedding layer is shared by both the generator and the discriminator.

Performer Encoder
Since its successful debut in Natural Language Processing (NLP) community, Transformer's attention mechanism has become a common recipe adopted across different domains of machine learning including speech processing [20], computer vision [1,9], and AIEd [23,2,11,24]. Compared to Recurrent Neural Network (RNN) family models, Transformer's attention mechanism has benefits of capturing longer-range dependencies and allowing parallel training, which enables the model to achieve better performance with less training time. However, despite these advantages, the time and memory complexities of computing the attention grow quadratically with respect to input sequence length, requiring demanding computing resources for training the model on long sequences. For instance, if L is input sequence length and d is dimension of query, key, and value vectors, Transformer's attention is computed as follows: where Q, K, V ∈ R L×d . The time and memory complexities for computing QK in the above equation are O(L 2 d) and O(L 2 ), respectively. Therefore, the cost for training Transformer becomes prohibitive with large L, preventing training the model even on a single GPU.
The problem of improving the efficiency of Transformer's attention mechanism is a common concern of machine learning community. Recent studies have proposed several methods to reduce the computing complexities lower than the quadratic degree with respect to input sequence length [19,30,17,27,5]. In this paper, we adopt Performer [5] since it uses reasonable memory and makes a better trade-off between speed and performance [28]. Performer approximates attention kernels through Fast Attention Via positive Orthogonal Random features (FAVOR+) approach. The R+part in FAVOR+ computes Q and K by applying random feature map φ : R d → R r + to each query vector q and key vector k in Q and K, respectively.
The FA-part leads to efficient attention mechanism computed as follows: where 1L ∈ R L is an all-ones vector of length L and diag is a diagonal matrix with the input as the diagonal. Unlike Transformer's attention mechanism which computes QK , applies the softmax, and multiplies the result and V , FA-VOR+ first computes K V , and multiplies the result and Q , reducing the quadratic time and memory complexities to linear ones with respect to input sequence length. Also, FA-VOR+ can model most attention kernels used in practice. If φ is a function given as below, FAVOR+ approximates the softmax attention kernel used in Transformer's attention mechanism.
Furthermore, the authors proposed a generalized attention kernel, an attempt to model a kernelizable attention mechanism beyond the commonly used softmax attention kernel, by defining φ as below.
where ReLU is the widely used Rectified Linear Unit function. In this paper, we use the generalized attention kernel instead of the softmax attention kernel since the former showed better performance in [5] and our preliminary experiments. Lastly, the O-part further reduces the variance of the estimator by making different random samples ws to be orthogonal to each other through the standard Gram-Schmidt orthogonalization process. For those who want to know more about theoretical details of FAVOR+, please refer [5].
With the efficient attention mechanism by FAVOR+, we propose the Performer encoder which is stacks of several identical reversible layers described in Figure 3. The reversible layer is based on Reversible Transformer [12,19] architecture to further improve memory efficiency in backpropagation. An input of the reversible layer x ∈ R L×d hidden is first chunked to x1, x2 ∈ R L×d hidden /2 . Then, scaled l2 normalization (ScaleNorm) [22] and FAVOR+-based multihead attention layer (MultiHeadAttn) are applied to x2, and the result is added to x1 to compute y1 ∈ R L×d hidden /2 . y1 = x1 + MultiHeadAttn(ScaleNorm(x2)).
After that, the scaled l2 normalization and point-wise feedforward layer (FeedForward) are applied to y1, and the result is added to x2, computing y2 ∈ R L×d hidden /2 . y2 = x2 + FeedForward(ScaleNorm(y1)). An output of the reversible layer y ∈ R L×d hidden is a concatenation of y1 and y2. We stack the reversible layer multiple times to allow the final model to sufficiently represent underlying data distribution.

Generator
The generator computes hidden representations [h G 1 , . . . , h G T ] by feeding the masked interaction sequence I M to a series of the interaction embedding layer (InterEmbedding), a point-wise feed-forward layer (GenFeedForward1), the Performer encoder (GenPerformerEncoder), and another pointwise feed-forward layer (GenFeedForward2): where I M E t , h G t ∈ R d emb and h GF t , h GP t ∈ R d gen hidden . Then, depending on whether the masked features are categorical or continuous variables, generator outputs are computed differently. If the masked features are categorical variables, the outputs are sampled from a probability distribution defined by the following softmax layer: If the masked features are continuous variables, the outputs are computed by the following sigmoid layer: . Similar to the case of categorical masked features, one can sample the outputs from a probability distribution defined by I M and parameters of the generator when the masked features are continuous variables. For instance, the outputs can be sampled from the Gaussian distribution where the mean and the variance are determined by I M and the generator's parameters. However, we make the outputs deterministic because sampling the outputs underperforms in our preliminary experiments when the masked features are continuous variables.

Discriminator
, and the sigmoid is applied to the last layer of the discriminator. After the pre-training, we slightly modify the discriminator by replacing the last layer with a layer having appropriate dimension for academic test performance prediction.

Training Objectives
The objective for pre-training is to minimize the following loss function: where GenLoss is the cross entropy (or mean squared error) loss function if the masked features are categorical (or continuous) variables, DisLoss is the binary cross entropy loss function, and 1 is the identity function. For ease of notation, we omit an index for each input sample in the above equation. If there are more than one masked features in each time-step (n > 1), the generator is trained under the multi-task leaning scheme. The objective for fine-tuning is to minimize the mean squared error loss between the model's predictions and score labels.

EXPERIMENTS 5.1 Dataset
The pre-training dataset consists of student interaction logs. The statistics of the dataset are summarized in Table 1. We exclude student interaction logs less than 15 in length to reduce noisy interactions from students using Santa just for trying out. Since most of the values of elapsed time, exp time, and inactive time are distributed in head areas as shown in Figure 4, we set their maximum to 300, 300, and 86400 seconds, respectively, and any values more than that are capped off to the maximum. We further normalize the values so that they are between 0 and 1 by dividing them by the maximum.
The fine-tuning dataset consists of test scores and student interaction logs before the test. Table 2 summarizes the statistics of the dataset. The number of score labels is far less than that of interactions in the pre-training dataset, which leads to the label-scarcity problem in academic test performance prediction. The minimum length of interactions before the test is 10 because Santa required students to solve at least 10 exercises before taking the test. The mean and median length of interactions before the test are longer than those in the pre-training dataset since students who decide to report their scores tend to be more serious about studying with Santa. As shown in Figure 4, most of the scores are distributed over the range from 700 to 900, and there are very few scores in the below 200 area. This has to do with the distribution of students using Santa whose initial and goal scores are usually in the range of 600 to 700 and higher than 800, respectively. Since the score labels are few in number, we perform 5-fold cross-validation by dividing the fine-tuning dataset into 5 splits and using 3/5, 1/5, and 1/5 of the dataset for training, validation, and test, respectively.

Training Details and Hyperparameters
We use the Mean Absolute Error (MAE) as the metric for academic test performance prediction. The list of hyperparameters and their values are described in Table 3. For each pre-training evaluation, the pre-trained model is fine-tuned and cross-validated on the validation set, which results in the same number of validation results as the number of pretraining evaluations. Then, we select the model with the best validation result and report an evaluation result of the model on the test set.

Effects of Generator's Pre-training Tasks
There are multiple interactive features to be masked in each token of the interaction sequence, which raises a question of how to construct a set of masked interactive features, and accordingly, which pre-training task for the generator is the most effective for academic test performance prediction. By default, all interactive features listed in Section 3 are taken as inputs for both the generator and discriminator. However, if response (or elapsed time) is masked, correctness (or timeliness) is excluded from the inputs and vice versa since there is an overlap of information that the features represent. For example, when both response and correctness are taken as inputs, and correctness is masked, the generator can predict the masked correctness by only looking at eid and response without considering other interactions, which leads to poor pre-training. The results are described in Table  4. The best result was obtained under the pre-training task of predicting response alone, which is slightly better than that of predicting correctness, and both response and correctness. Predicting correctness of student response is an important task in AIEd as can be seen from the large volume of studies about Knowledge Tracing. Also, [3] empirically showed that student response correctness is the most pedagogical interactive feature for academic test performance prediction. However, rather than pre-training a model to predict whether a student correctly responded to a given exercise, the pre-training task of predicting student response itself injects more fine-grained information into the model, which leads to the more effective pre-training for academic test performance prediction. Interestingly, the underperformed results were obtained when predicting elapsed time or timeliness in pre-training despite the benefits their information bring to several AIEd tasks [10,33,24]. We hypothesize that elapsed time and timeliness may introduce irrelevant noises and thus guide the model towards a direction inappropriate for academic test performance prediction. In the case of exp time and inactive time, we observed that the generator failed to learn to predict their values when only given the interactive features listed in Section 3, which leads to unstable pre-training. From these observations, in the following subsections, we conduct experimental studies based on the pre-training task of predicting response alone.

DPA vs. Baseline Methods
We compare DPA with the following pre-training methods: • No pre-training: We train the fine-tuning models only on the fine-tuning dataset.
• Autoencoding: Autoencoding (AE) is a generative pretraining method widely used across different domains of machine learning including AIEd [13,8]. Given an unmasked interaction sequence, AE pre-trains a model to reconstruct the input interaction sequence.
• Assessment Modeling: Assessment Modeling (AM) [3] is the previous state-of-the-art generative pre-training method for academic test performance prediction. In AM, a model takes a masked interaction sequence as an input and is pre-trained to predict masked features. AM is exactly the same as fine-tuning the pre-trained generator in DPA.
Also, we investigate whether DPA is effective with the following different fine-tuning models: • MLP: Multi-Layer Perceptron (MLP) is stacks of simple fully-connected layers. Given an interaction sequence, interaction embedding vectors of all time-steps are summed together to compute a fixed-dimensional vector which is fed to a series of the fully-connected layers.
• BiLSTM: Bi-directional Long Short-Term Memory (BiL-STM) is a model widely used for time series data prediction tasks. The global max pooling layer is applied on top of the BiLSTM layer to obtain a fixeddimensional intermediate representation from an input sequence of varying length.
• Transformer Encoder: Transformer Encoder is a series of several identical layers composed of a multi-head self-attention layer with the softmax attention kernel and a point-wise feed-forward layer. We set the Transformer encoder's attention window size to 512 due to the out of GPU memory occuring when training the Transformer encoder of 1024 attention window size on our single GPU machine.
As described in Table 5, transferring the pre-trained knowledge brings better results in most cases, and the best result is obtained from DPA. Especially, when the Performer encoder, the best performing fine-tuning model, is used as the fine-tuning model, DPA reduces MAE by 4.05%, 21.84%, and 38.05% compared to AM, AE, and No pre-training, respectively. Among the baseline pre-training methods excluding No pre-training, the worst result is obtained from AE beacuse the pre-training task of AE is much easier than  that of AM and DPA. We observed that the loss curve of AE converged to near zero within the first pre-training evaluation.

Robustness to Increased Label-scarcity
Since the motivation behind our proposal of DPA is the label-scarcity problem, we investigate how MAE changes at varying degrees of label-scarcity. Figure 5 and Table 6 describe the results when using 1/2, 1/4, and 1/8 of the total number of fine-tuning training samples. In all degrees of label-scarcity, DPA consistently outperforms AM. Also, DPA fine-tuned on 1/2, 1/4, and 1/8 of the dataset outperforms AM fine-tuned on the entire dataset, 1/2, and 1/4 of the dataset, respectively, which shows that DPA is more robust to label-scarcity than AM. Compared with No pre-training, the gap between No pre-training and the other two pre-training methods increases as the number of labels becomes scarce. Furthermore, the other two pre-training methods fine-tuned on 1/8 of the dataset outperform No pre-training fine-tuned on the entire dataset. Figure 4 described that score labels are mainly distributed over the specific ranges. We investigate how this biased distribution of score labels affects model predictions. The results are described in Figure 6. As expected, DPA severely underperforms when the score labels are below 200. Although these are natural results from a machine learning perspective, this is a serious problem from a perspective of educational service because students whose scores are lower than 200 are inaccurately diagnosed their academic status. It is also against the equity of education since not all students can receive the same level of educational service. There may be various research directions to solve this problem, such as generating pseudo labels, measuring prediction uncertainties, or even collecting more score labels. However, we don't go deeper into it any further and leave it as a future work.

Analysis of DPA
The previous experimental results showed that DPA makes better predictions and is more robust to label-scarcity than AM. However, where the gains from DPA are coming from is not obvious. We investigate what makes DPA outperform AM by comparing DPA and AM with the following set of ablation pre-training methods: • DPA60%: DPA60% is the same as DPA except the discriminator loss tokens only come from the masked interactions. Since we set the masked interaction ratio to 60%, DPA60% pre-trains the discriminator with 40% fewer loss tokens than DPA.
• RAM: Unlike DPA where there are the generator and the discriminator, in Replaced Assessment Modeling (RAM), there are two generators of different sizes, a small and large generator. The small generator is pretrained in the same way as AM. The large generator takes a replaced interaction sequence which is generated by replacing the masked features with the small generator's outputs as an input, and is pre-trained to predict the masked features. After the pre-training, we throw away the small generator and fine-tune the large generator. The sizes of the small and large generator are the same as those of DPA's generator and discriminator, respectively.  • AAM: All-tokens Assessment Modeling (AAM) is the same as RAM except the large generator is pre-trained to predict all features rather than just predicting the masked features. Figure 7 depicts graphical description of each pre-training method. The results are described in Table 7. In the following subsections, we analyze the results in aspects of pretrain/fine-tune discrepancy due to the mask token, discriminative vs. generative pre-training, and sample efficiency.

Pre-train/Fine-tune Discrepancy Due to Mask Token
During pre-training, the generator in AM sees the mask token which does not appear in fine-tuning, leading to pretrain/fine-tune discrepancy. [31,6] also raised the same issue found in the MLM pre-training method proposed in [7]. Since the discriminator in DPA does not see the mask token both in pre-training and fine-tuning, DPA does not suffer from the pre-train/fine-tune discrepancy, making it necessary to examine how much gain is obtainable from this issue. Comparing RAM and AM, removing the pre-train/fine-tune discrepancy slightly reduces MAE by 0.28%.

Discriminative vs. Generative Pre-training
The discriminator in DPA and the generator in AM are pretrained with different objectives. For each token to be pre-dicted, the discriminator pre-training loss comes from a discrimination error between the discriminator's output and an originality of the token. On the other hand, the generator is pre-trained by generation error between the generator's output and an identity of a token to be predicted. Comparing DPA with AAM, and DPA60 with RAM, the discriminative pre-training objective reduces MAE by 0.82% and 1.46%, respectively, over the generative pre-training objective.

Sample Efficiency
The loss tokens for pre-training the discriminator in DPA come from all interactions in the input interaction sequence. However, the generator in AM is pre-trained with loss tokens only coming from the masked interactions. Considering the masked interaction ratio is 60%, DPA is more sample efficient than AM since the discriminator is pre-trained with 40% more loss tokens than the generator. Comparing DPA with DPA60, and AAM with RAM, sample efficient pre-training brings reduction in MAE by 2.35% and 2.89%, respectively, showing that the sample efficiency contributed most to the improvement from AM to DPA.

CONCLUSION
In this paper, we proposed DPA, a transfer learning framework with discriminative pre-training tasks for academic performance prediction. Our experimental results showed the effectiveness of DPA for the label-scarce academic performance prediction task over the previous state-of-the-art generative pre-training method. Avenues of future research include investigating more effective pre-training tasks for academic performance prediction and pre-train/fine-tune relations in AIEd.


Virtual Training Environment for Gas Operatives: System Usability and Sense of Presence Evaluation

Training of gas operators in real-life settings often has associated risks to health and property. The use of a virtual environment to train gas operators has the potential to offer risk-free training. This study tests the usability of a virtual environment specifically designed to teach new gas operatives in near real-life scenarios. Thirty-two participants tested the virtual environment and performed different tasks required to complete gas safety checks. We used SUS (System Usability Scale) and sense of presence questionnaires to collect data from these participants. The SUS analysis indicated that most participants belonging to a different gender, age, and virtual reality experience groups were comfortable in the VR training environment. The sense of presence data analysis also confirmed similar results as all sense of presence factors scored high regardless of the demographics characteristics of the participants. However, there is still a need to add different scenarios to make the virtual environment into a comprehensive training course.

INTRODUCTION
The recent statistics show that carbon monoxide (CO) poisoning has resulted in thousands of incidents. These incidents have caused more than 697 deaths inside the UK in the last two decades (Safety, 2018). The biggest reason for CO emissions is improper installation of gas ovens, boilers, and cookers. Lack of maintenance of such appliances also contributes to such incidents.
Usually, in gas leak incidents, the gas operatives are called in to fix such situations. Therefore, proper training of gas operatives is imperative to avoid any causalities. However, training new gas operatives in real-world settings has its associated risks. Developing a full-scale gas leak training scenario will not only come at a high cost but will have associated risks to life and property as well. Considering these factors, the researchers have previously advised incorporating virtual reality (VR) as a potential solution to reduce these risks to life and property (Asghar, Egaji, Dando, Griffiths, & Jenkins, 2019). Additional research has also shown that incorporating VR in traditional training can provide a more interactive environment for the trainees and positively impact their skill set and knowledge retention (Fredricks, Blumenfeld, & Paris, 2004).
In literature, we can find many examples of virtual environments to help train people in risk-free environments. Examples include VR training for firefighters in fire evacuation (Cha, Han, Lee, & Choi, 2012). A simulator for tunnel fire evacuation (Ronchi et al., 2015). VR-based emergency rescue training system for railway cranes operators (Xu et al., 2018). VR application for training children to cross railway crossings safely (Dando, Asghar, Egaji, Griffiths, & Gilchrist, 2018). Another VR application to train miners in mine incidents (Kizil & Joy, 2001). These examples motivated us to develop a virtual environment for safety training of the gas operatives explained our previous work (Asghar et al., 2019).
The current paper study aims at analysing the usability of a virtual training environment for the gas operatives. Two data collection instruments called system usability scale (SUS) and sense of presence questionnaires are used to test the usability and level of immersion of the gas operatives while trained through the virtual training environment. In total, 32 participated in the virtual environment testing across multiple sessions. Most of the participants supported the idea of the virtual environment and appreciated the skills and knowledge this system offered to them. 140 The rest of the paper is divided into five sections. A summary of the virtual training environment is presented in the second section. The research process used for this study is summarised in the third section. Research results and discussions are shown in the fourth section. Finally, conclusions from this study and future work recommendations are offered in the fifth section.

VIRTUAL TRAINING ENVIRONMENT
The virtual training environment is based on two scenarios. The first scenario presents a typical residential home in virtual settings. The gas operatives are expected to explore any gas or CO emitting appliances in the virtual home and select these as potential hazardous appliance for gas leaks. The gas operatives can use teleporting devices to enter the house. Teleporting helps them open or close the doors, select appliances, take their readings, and mark any appliance deemed dangerous. The gas operatives can select multiple appliances within one scenario. At the end of this scenario, the gas operatives can see their results in the form of a summary. The gas operatives can repeat this scenario many times. The first scenario is presented in figure 1.

Figure 1. First Scenario with Potential Actions
The second scenario requires the gas operatives to follow a storyboard of a gas leak situation, which can occur in real-home settings. In this scenario, the gas operative will have a time window of 40 seconds to follow the systematic approach from hazards' identification, inspection, and solution. The second scenario showing potential actions that a gas operative can take in a particular situation is presented in figure 2.
The gas operative will read a particular appliance, mark that appliance as hazardous or disconnect that appliance. If the gas operative cannot find the gas leak in 40 seconds or commits any mistakes in this process, there will be an explosion in the virtual training environment. This will add to the realness of dangers involved in such situations in real-life settings. Both scenarios are expected to aid gas operatives and be an excellent addition to their existing training courses.

Figure 2. Second Scenario with Potential Solutions
The current virtual training environment is kept limited to two scenarios only, as previous research has shown that excessive VR exposure can result in dizziness and sickness for some users (Viirre & Ellisman, 2003).

RESEARCH PROCESS
This study followed a three-tier research strategy to accomplish the main research aim as presented in figure 3. Tier one focused on the design of the virtual training environment; tier two concentrate on testing and data collection from the gas operatives, and tier three concentrate on the analysis of the collected data. The ethical approval for the study was obtained from the Faculty of Computing, Engineering and Science, University of South Wales.

RESULTS AND DISCUSSION
In total, 32 participants were trained in the virtual training environment. The training process involved each participant going through both scenarios in the virtual world and filling out both questionnaires at the end of their experience. Table 1 summarises the demographic information for the participants. We have followed the standard SUS data analysis process for data analysis and result interpretation (Harrati, Bouchrika, Tari, & Ladjailia, 2016). For the sense of presence questionnaire data, we used the Shapiro-Wilk and the Mann-Whitney U tests.

SUS Data Analysis
The average historical score for SUS is 68. A SUS score greater than 68 means that the system usability is good, and a lower than average score means there is some problem with the usability of that system (Bangor, Kortum, & Miller, 2008). For the current study, 29 participants have SUS scores of above 68, which means most participants were satisfied with the usability of the virtual training environment. Only three participants with scores of less than 68 were not fully satisfied with the usability of the virtual training environment.
Another dimension in SUS analysis is known as grading levels of SUS scores. In such research, a score equal to or greater than 80.3 means grade A, 68 and above is grade C, and equal to or less than 51 means grade F (Sauro, 2011). The average SUS score of all 32 participants who used the VR training environment is 85.31. This overall survey score corresponds to grade A, which means that the participants did enjoy using the virtual training environment, and they will recommend it to others. Table 2 shows SUS scores based on the demographics of the participants. Interestingly, all groups scored grades A, which confirms the quality of the virtual training environment.

Sense of Presence Data Analysis
For the sense of presence data collection and analysis, we have considered the work of Witmer & Singer, as they are believed to be the pioneers in this domain. Their work revolved around the theories of involvements and immersion. Based on their theoretical and empirical research work Witmer & Singer determined multiple factors, including control, sensory, distraction, and realism, contribute to the sense of presence and immersion (Witmer & Singer, 1998).

Factors Impact on Sense of Presence
The sense of presence was measured on a 7-point Likert scale. The basic descriptive statistics show mean values for control (6.16), sensory (5.97), distraction (4.06), and realism (5.44). The higher values of control, sensory, and realism indicate that sense of presence in the virtual training environment was high amongst the participants. The low value for distraction suggests a minor degree of interference by objects and actions in the virtual training environment, and participants could concentrate on the tasks. In summary, all four factors in this study contribute to a high sense of presence in the virtual environment.

Demographics Impact Sense of Presence
For the second part of the analysis, we wanted to test the impact of demographic characteristics on the presence factors. The Shapiro-Wilk test showed that all demographic characteristics have significant values (p< 0.001), which means the data set in this study was not normally distributed. Therefore, we have to use some non-parametric tests for testing differences among different data groups. As all demographic characteristics, including gender, age group, and previous VR experience, have only two possible values, we used the Mann-Whitney U test for further analysis (Nachar, 2008).
The results of the Mann-Whitney U test are summarised in table 3. These results show no significant values of (p< 0.05) for any demographic characteristics of the participants. This indicates that none of these characteristics has a significant impact on the system's usability under test. Therefore, we can conclude that all participants, regardless of their gender, age, and VR experience, enjoyed the virtual training environment as the factor mean scores indicated in the last subsection.
From the data and results for both questionnaires, we can conclude that system usability and sense of presence are highly appreciated by the gas operatives who got the chance to have hands-on experience in the virtual training environment. Furthermore, the usability results are slightly better than the previous study (Asghar et al., 2019). This indicates that involving end-users in the testing process and modifying the system based on their feedback results in improved system usability and sense of presence. The study limitation includes that this system is still undergoing further development based on user feedback. Furthermore, as only 32 people participated in the testing, the results cannot be generalised and need further testing.

CONCLUSION AND FUTURE WORK
A virtual training environment for the gas operatives is evaluated in this paper. This virtual training environment can be integrated into existing modules of gas operatives training that can help in improving the skill set and decision-making capabilities among the gas operatives. The potential advantages of using a virtual training environment are its ability to construct close to reallife situations at less cost, risk-free training of new gas operatives, and execution of the same problem multiple times.
The virtual training environment is tested with 32 participants with the help of the SUS and sense of presence questionnaires. The average SUS score for all participants was 85.31, equivalent to a Grade A in the SUS analysis. There are no significant differences between different age groups, gender, and participants with or without previous VR experience. The sense of presence data analysis supports the results from the SUS scores reflecting high system usability and close to a real-life virtual training environment.
Future work will address some of the limitations of the current paper by testing the virtual training environment with a much larger sample size to generalise the results. In addition, more VR scenarios can add extensive details and interaction opportunities between the gas operatives and their trainers by having trainee and instructor views (multiplayer).