## 6 Conclusions and Discussion
\begin{table}
\begin{tabular}{c r r r r r r r r r r r r r} \hline \hline \(\theta\) & \(\alpha\) & \(\beta_{1}\) & \(\beta_{2}\) & \(\beta_{3}\) & \(\beta_{4}\) & \(\sigma_{A}^{2}\) & \(\sigma_{B11}^{2}\) & \(\sigma_{B22}^{2}\) & \(\sigma_{B33}^{2}\) & \(\sigma_{B44}^{2}\) & \(\sigma_{e}^{2}\) & \(\pi\) & \(\lambda\) \\
**True** & \(\mathbf{0.42}\) & \(\mathbf{10.52}\) & \(\mathbf{3.89}\) & \(\mathbf{2.72}\) & \(\mathbf{-0.10}\) & \(\mathbf{0.34}\) & \(\mathbf{4.23}\) & \(\mathbf{0.80}\) & \(\mathbf{0.70}\) & \(\mathbf{0.50}\) & \(\mathbf{0.11}\) & \(\mathbf{0.32}\) & \(\mathbf{1.09}\) \\ \hline \multicolumn{11}{c}{} \\ Est & \(0.38\) & \(10.39\) & \(3.87\) & \(2.68\) & \(-0.47\) & \(0.37\) & \(5.82\) & \(0.90\) & \(1.42\) & \(0.50\) & \(0.12\) & \(0.32\) & \(1.07\) \\ MSE & \(0.33\) & \(0.03\) & \(0.04\) & \(0.05\) & \(0.39\) & \(0.39\) & \(0.44\) & \(0.23\) & \(1.14\) & \(0.25\) & \(0.56\) & \(0.33\) & \(0.18\) \\ Bias & \(-0.04\) & \(-0.13\) & \(-0.02\) & \(-0.04\) & \(-0.30\) & \(0.03\) & \(1.59\) & \(0.01\) & \(0.72\) & \(0.00\) & \(0.07\) & \(-0.00\) & \(-0.02\) \\ \multicolumn{11}{c}{} \\ \multicolumn{11}{c}{} \\ Est & \(0.42\) & \(10.23\) & \(3.78\) & \(2.23\) & \(-0.57\) & \(0.97\) & \(8.16\) & \(1.08\) & \(1.06\) & \(1.00\) & \(0.15\) & \(0.29\) & \(1.15\) \\ MSE & \(0.45\) & \(0.04\) & \(0.05\) & \(0.60\) & \(3.50\) & \(1.93\) & \(0.94\) & \(0.32\) & \(0.83\) & \(1.25\) & \(0.57\) & \(0.34\) & \(0.21\) \\ Bias & \(0.00\) & \(-0.29\) & \(-0.11\) & \(-1.79\) & \(-1.80\) & \(0.63\) & \(3.93\) & \(0.25\) & \(0.42\) & \(0.59\) & \(0.09\) & \(-0.03\) & \(0.06\) \\ \hline \hline \end{tabular}
\end{table} 
Table 3: Simulation results on the performance of StEM algorithm when the post-ART trajectory is modeled with two-compartment model (M2). more efficient variants of the algorithm, such as independent sample approaches [32], to accommodate these extensions. Several methodological developments offer promising directions. One is the integration of machine learning (ML) techniques--for example, to model post-ART viral dynamics or improve scalability. While ML algorithms may enhance flexibility, they often lack the interpretability and inferential tools required for surveillance-oriented estimation tasks. Hybrid approaches combining ML with structured statistical models could offer the best of both worlds and merit future study. Another important extension involves incorporating subject-level covariates (e.g., age, gender, or transmission risk) into the fixed effects components of the model. This would allow parameters such as \(\alpha_{1}\), \(\pi\), \(\lambda\), \(\beta_{1}\), \(\beta_{2}\), and \(\beta_{3}\) to vary across individuals. However, this generalization would require replacing the closed-form M-step with iterative procedures such as Newton-Raphson, significantly increasing computational demands--particularly in the presence of censoring and complex random effects. We are actively pursuing this extension using parallel computing strategies to support large-scale applications. It will also be important to assess the robustness of our method under misspecified change point distributions. While we focused on the ZIE distribution here, the algorithm can, in principle, accommodate zero-inflated log-Gaussian, gamma, or Weibull alternatives, as mentioned above. Incorporating these would, however, increase the complexity of Metropolis-Hastings updates and reduce sampling efficiency. Simulation-based robustness assessments under alternative distributions are planned as the next step. Conceptualization, S.L.B., D.N. and H.Z. ; methodology, H.Z. ; preliminary data preparation, M.R. ; data curation, D.B.H. and U.R.F. ; writing--original draft preparation, review and editing, L.W. and H.Z. All authors have read and agreed to the published version of the manuscript. This work was partially supported by NIH grant R21 AI147933. It was also partially supported by the High-Performance Computing Center at the University of Kentucky.