I am trying to make a discrete SEIR model with the I category partitioned into two; Iasymptomatic and Isymptomatic, thus it is a SEIIR model. I solve it numerically using Euler method (in python). I've read a lot about the theory on discrete SIR, but it seems like I can't apply the theory in my model and my question is why...
My model is the following difference equations: \begin{align} ∆S&=−βS(δI_a+I_s)\\ ∆E&=βS(δI_a+I_s)−ψE\\ ∆I_a&=(1−η)ψE−γI_a\\ ∆I_s&=ηψE−γI_s\\ ∆R&=γ(I_a+I_s) \end{align} where $β=$ transmission rate, $γ=$ recovery rate, $η =$ share of population symptomatic, $ψ =$ exposed rate, $(1-δ) =$ indicates the percentage decrease to transmit infection for asymptomatic individuals.
I've read that the maximum of $I(t)$ in a discrete SIR-model happens when $S(t) = (1/R_0)S(0)$ where $R_0$ is the basic reproduction number. However, this happens not to be the case for me when I use $I(t) = I_a(t)+I_s(t)$. (I calculated $R_0$ by using next generation metrix and I am pretty sure my result here is correct)..
My questions are: Can I use the theory for SIR to my SEIIR model? How do I find the exact solution to this system? I have an idea that I want to find the error in my numerical solution, but that means having the exact solution at hand, right?
Any help is much appreciated ! Thanks in advance