Fitting parameters together with initial conditions of an ODE to real data

Yes

Is there any better way of setting?

You can also use remake and set parameters and initial conditions accordingly like prob = remake(prob; p = p, u0 = X0)

How about controlling the bound of parameters based on other parameters? Is this truncated(Normal(0,Λ*(ηA+ηS)/(μ*μp)); lower=0, upper=min(Λ*(ηA+ηS)/(μ*μp),65188)) too dumb?

I think that’s the way to do it