## 2 Materials and Methods
### Simulation Framework  
The simulation of the bioequivalence in the context of 2 \(\times\) 2 crossover design consisted of generating subjects for the reference (i.e., R) and the test (i.e., T) groups for both periods, termed as \"original\" [6]. Initially, a total of N\({}_{\text{R}}\) subjects was generated for the first period, corresponding to the R product, from a random process with mean \(\upmu_{\text{R1}}\) and standard deviation \(\sigma_{\text{R1}}\). For the T drug, N\({}_{\text{T}}\) subjects, equal to N\({}_{\text{R}}\), were generated, following the same process with mean \(\upmu_{\text{T1}}\) and standard deviation \(\sigma_{\text{T1}}\). Both groups had equal coefficient of variation (CV), referring to the between-subject variability. For the second period, each of the total N = N\({}_{\text{R}}\) + N\({}_{\text{T}}\) subjects was multiplied with a stochastic term, which was generated through a random process, to incorporate the within-subject variability of the subjects. The original data were randomly subsampled with varying proportions, referred to as \"subsampled\". These subsampled datasets were used to train the VAE model and, in a next step, to synthetize the new virtual subjects termed as \"generated\". Finally, the \"generated\" data from both groups (T and R) and study periods (i.e., periods I and II of the 2 \(\times\) 2 crossover design) were transferred for statistical analysis, using the official statistical framework imposed by regulatory authorities for bioequivalence [4; 5]. This statistical analysis refers to the following steps: ln-transformation of the variables, application of a linear model (ANOVA), calculation of the residual error, using this error term to construct a 90% confidence interval for the mean difference between T and R in the ln-domain, and, lastly, checking if BE is declared [4; 5; 6]. This entire procedure was repeated 500 times and the percentage of bioequivalence acceptance was measured at the end [6]. The above-mentioned route of analysis can be outlined in the following steps:  
1. N individuals are randomly generated for both groups in the case of the first period: 1. N\({}_{\text{T}}\) individuals for the T group for the first period, with mean \(\upmu_{\text{T1}}\) and standard deviation \(\sigma_{\text{T1}}\); 2. N\({}_{\text{R}}\) individuals for the R group for the first period, with mean \(\upmu_{\text{R1}}\) and standard deviation \(\sigma_{\text{R1}}\); 3. The T and R groups were set to have equal CVs; 4. The sample sizes of T and R groups were assumed to be equal: N\({}_{\text{T}}\) = N\({}_{\text{R}}\); 5. Thus, the sample size of the study is: N = N\({}_{\text{T}}\) + N\({}_{\text{R}}\). 2. The N individuals from the first period are multiplied, with a randomly generated \"stochastic term\", with mean \(\upmu_{\text{ST}}\) and standard deviation \(\sigma_{\text{ST}}\). The stochastic term coefficient of variation (CVw) represents the \"within subject variability\" for each simulated volunteer between the period I and II of the crossover study [6]. It should be noted that the CVw is different to the between-subject variability (i.e., CV) discussed in step \"Ic\". \begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline \multicolumn{2}{c}{**Activation Function**} & \multicolumn{2}{c}{**Weights of Loss Function**} & \multicolumn{2}{c}{**Number of Hidden Layers**} & \multicolumn{2}{c}{**Number of Neurons in Hidden Layers**} \\ \hline
**Hidden Layers** & **Output** & **Kullback– Layer** & **Reconstruction** & **Encoder** & **Decoder** & **Encoder** & **Decoder** \\ \hline Softplus & Linear & 1 & 1 & 3 & 3 & 64-32-16 & 16-32-64 \\ \hline  & & 2 & 2 & & & & \\ \hline  & & \(\dots\) & \(\dots\) & & & & \\ \hline  & 9 & 9 & & & & & \\ \hline  & 10 & 10 & & & & & \\ \hline \hline \end{tabular}  
* The entire computational work and “runs”, outlined in Table 1, were performed with TensorFlow 2.10.0 and Python version 3.7 and executed within a “Jupyter notebook” environment.  
\end{table}