The aim of this paper is to model the dynamics of the human papillomavirus (HPV) in cervical epithelial cells. We developed a mathematical model of the epithelial cellular dynamics of the stratified epithelium of three (basale, intermedium, and corneum) stratums that is based on three ordinary differential equations. We determine the biological condition for the existence of the epithelial cell homeostasis equilibrium, and we obtain the necessary and sufficient conditions for its global stability using the method of Lyapunov functions and a theorem on limiting systems. We have also developed a mathematical model based on seven ordinary differential equations that describes the dynamics of HPV infection. We calculated the basic reproductive number (
R
0
) of the infection using the next-generation operator method. We determine the existence and the local stability of the equilibrium point of the cellular homeostasis of the epithelium. We then give a sufficient condition for the global asymptotic stability of the epithelial cell homeostasis equilibrium using the Lyapunov function method. We proved that this equilibrium point is nonhyperbolic when
R
0
=
1
and that in this case, the system presents a forward bifurcation, which shows the existence of an infected equilibrium point when
R
0
>
1
. We also study the solutions numerically (i.e., viral kinetic in silico) when
R
0
>
1
. Finally, local sensitivity index was calculated to assess the influence of different parameters on basic reproductive number. Our model reproduces the transient, acute, latent, and chronic infections that have been reported in studies of the natural history of HPV.