How to perform parameter estimation for ODEs, including simulating steady steady state and using multiple experimental conditions in Julia

Question

I am trying to figure out how to do parameter estimation of a system of ODEs for multiple experimental conditions, with unknown initial values. I know the relationships between the different quantities of the different variables, but not the exact values. Lets say we have a reaction going between R <--> Rp, influenced by some parameters. We know R + Rp will sum to 1, but we do not know if R = 1, Rp = 0, or R = 0.7, Rp = 0.3. Furthermore, the system can be given some known stimulus (S). For these different stimuli, experimental data is available, for one variable (Rp). I would like to be able to take all of this into account and estimate the parameters that gives simulations closest to the experimental data. The current approach implemented in MATLAB/Python (is to simulate what the steady state would be, given some arbitrary values (based on the constraint mentioned above), with no stimuli. Then, the values from that steady state simulation is used as initial values for the next simulations. In the end, all simulations should be compared to experimental data (e.g. with an L2-term), summed for all experiments.

Tl, dr; I want to be able to:

1. Simulate steady state
2. Use the steady state as initial values for simulation with some input (e.g S=1)
3. Use steady state for simulation with another input (e.g S=2)
4. Calculate loss (summed) for both simulations with corresponding datasets.
5. Optimize the model parameters In my current code (see below) I can find a good parameter estimation for one dataset, if I have a pre-defined guess of the steady state. When I try use the SteadyStateProblem() function, I do not get a reasonable estimation of the steady state values. How can I extend my script to also be able to use the simulated steady state, and optimize for multiple experimental conditions?

Thanks in advance!

## Define packages
using Plots
using DiffEqParamEstim
using Optim
using DifferentialEquations

## Setup data (for one condition), t = time, µ = mean values. SEM = data uncertaintes
t = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0]
µ = [-0.006511181, 0.0769646, 0.0601607, 0.0509726, 0.0376803, 0.0279487, 0.0349815, 0.0282609, 0.02870767, 0.0190136, 0.0218868] 
SEM = [0.0027889401688424983, 0.0035224393858044835, 0.0045910703207422115, 0.005134166212086762, 0.003473086632473579, 0.0036044191797477347, 0.0038831248191515955, 0.001998941711395196, 0.004765052514110766, 0.0036907553485251046, 0.001812117336156795]

## Define model
function M1!(du,u,p,t,input =[1]) # Represents M1. Note that some parameter orders have switched. 
    R, Rp, RS, RSp = u
    k1, k2, k3, k4, kfeed = p
    S = input[1]

    r1 = R*S*k1 
    r2 = Rp*k2
    r3 = Rp*RSp*kfeed 
    r4 = RSp*k3
    r5 = RS*Rp*k4

    du[1] = dR = r3+r2-r1 
    du[2] = dRp = r1-r2-r3 
    du[3] = dRS = r4-r5 
    du[4] = dRSp = r5-r4
end

## Enable variable input via enclosure
input = 0
f = (du,u,p,t) -> M1!(du,u,p,t,input) # Closure to enable different inputs

## Simulate steady state for initial parameterset p0
x0_good = [1.0 0 1 0] # initial values
x0_naive= [0.5 0.5 0.5 0.5] 
p0=[1 0.0001 1 0.01 1000000] #start guess

stead_state = solve(SteadyStateProblem(f,x0_naive, p0))
println("Solved SteadyStateProblem: $(stead_state)")
simulated_steady_state = solve(ODEProblem(f,x0_naive, (0.0, 100000), p0))
println("Simulated steady state: $(simulated_steady_state.u[end])")
println("\"Known\" initial values: $(x0_good)")

## Define ODE problem, setup objective function
input = 1 #Set the input S = 1
prob=ODEProblem(f, x0_good,  (minimum(t), maximum(t)), p0) 
loss_function = L2Loss(t, µ, data_weight=SEM)
obj = build_loss_objective(prob,Tsit5(), loss_function, saveat=t,save_idxs=[2], maxiters=10000)

## Optimize the model parameters (for one input condition, S=1)
lb = (ones(1, length(p0)) * 1e-7)
ub = (ones(1, length(p0)) * 1e7)
results = optimize(obj,lb,ub,p0) # Do the optimization. 
popt = results.minimizer

## Plot results
prob = ODEProblem(f, x0_good, (minimum(t), maximum(t)), popt) # Same as the prob above, but with a new parameter set. 
sol = solve(prob, saveat=t, save_idxs=[2]) 
solHR = solve(prob, save_idxs=[2]) 

plotly() # Select plotting backend
plot(solHR, label="sim, highres", linewidth=2) # Plot simulation with high time resolution
scatter!(sol, label="sim at experimental time points") # Plot simulation at experimental timepoints
scatter!(t,µ,yerror=SEM,label="data") # scatter data