Variational autoencoders are capable of replicating sequence statistics and capturing evolutionary trends. Before testing our hypothesis, we embarked on selecting the best model architecture based on the implementation provided by Ding et al. . This involved optimizing the encoder, decoder, and training procedure to minimize the difference between the generated and input sequence distributions (generative capacity) while also preserving the relationship between phylogeny and the latent space (geometric properties). In order to evaluate the model's generative capacity, we implemented several tests. The first test examined how well our model reproduces the statistics of the input dataset on the output. To this end, we compared the first and second-order statistics of 3,000 randomly sampled MSA input sequences with those generated by our VAE model (Fig. ). The first-order statistics compare the frequencies of amino acids in every position. The VAE model successfully reproduced them, also maintaining the frequencies of gaps, indicating that it did not introduce extra deletions. Although first-order statistics often provide strong signals for belonging to a given protein family, they cannot secure the generation of viable proteins, as the interaction between distant positions can be critical . The second-order amino acid distribution can serve as an important metric to evaluate the similarity of distant relations in the generated sequences. We employed secondorder covariances to assess the ability of the model to reproduce this distribution, as suggested by prior works . The VAE model reproduced the empirical second-order statistics of natural sequences (Fig. ). We integrated query reconstruction accuracy as an additional metric in our analysis to ensure that the model was capable of reconstructing the query sequence with minimal mutations. Notably, the final Model 1 demonstrated the ability to reconstruct the query with as few as 6 unknown mutations while trained on HLD I-IV dataset. In our second test, we examined the statistical profile of our model by measuring the model average reconstruction accuracy of sequences from various control sets. The average reconstruction accuracy of the sequence was approximated as an average reconstructed sequence identity for 5,000 samples around the original sequence coordinates of its latent space embedding. For the negative control, we generated a random set of sequences with matching frequencies of amino acids in every position to those in the input MSA. The training subset consisted of the same number of sequences sampled randomly from the training dataset. For the positive control, we evaluated the average reconstruction accuracy of 5% MSA sequences removed from the input before training. We ended up with 612 sequences for the negative control, positive control, and training subsets. The test showed that our model can distinguish random sequences from those in the MSA with the cutoff value of the reconstructed accuracy of 50% (Fig. ).