Hello All,
I have a coupled set of differential equations with sometimedependent parameters that I would like to optimize.
In the code at the bottom of this post, the k_i terms are fixed kinetic constants, the c[index] terms are concentrations of chemical species, and the s[index] terms are the source terms which are being solved for. S and c are timedependent. I also have observations at discrete timesteps (0, 1, 2, 3, …, n) of each species. My questions are the following:

For solving the forward problem, The parameters (sources) have values at fixed, in this case discrete, points in time. However, the ODE solver will take time steps between these discrete values (e.g., t = 1.2 versus 1 or 2). So far, I perform an interpolation for the source terms in the time domain by interpolating the parameters in time. Please refer to the function interpSources in the code. Is there a better way to get parameters that are continuous in time e.g., passing atimedependent function into ODEProblem instead of an array of parameters?

Also, do you have any recommendations on how to proceed with optimizing thesetimedependent parameters (sourceterms) constrained by observational timeseries? I tried touring_inference from DiffEqBayes.jl, but that doesn’t seem to work due to the interpolation problem mentioned in Question 1. The examples online are for fixed parameters (e.g., the Lawrence Equations or predatorprey model).
The simplified code is below. Please let me know if you need more clarification, because I know this is a rather convoluted question. Thanks in advance!
using DifferentialEquations, DiffEqBayes, ParameterizedFunctions
using Distributions, Interpolations
end_time = 25.0;
time = collect(1.0:end_time);
num_steps = length(time);
tspan = (1.0 ,end_time);
# define some constants
const k1, k2, k3 = 7.89, 464.8, 3600*24;
function interpSources(s, t)
num_time, num_species = size(s)
species_grid = [i for i =1:num_species];
grid = (time, species_grid);
intp = interpolate(grid, s, Gridded(Linear()));
s_out = intp(t, species_grid);
return s_out
end #function interpSources
function BoxModel(dc ,c ,p ,t)
s = interpSources(p ,t) # interpolate the source terms (parameters) in time
dc[1] = s[1]  k1 * c[1] * c[3];
dc[2] = s[2] + k1 * c[1] * c[3]  k2 * c[2] * c[3];
dc[3] = s[3]  k1 * c[1] * c[3]  k2 * c[2] * c[3]  k3 * c[3];
end # function BoxModel
# define initial conditions
c_0 = [1700.0 100.0 3.71e5];
# defining the sourceterms (parameters) where row i is time i and col j is species j
sources = [1.05 * ones(num_steps, 1) 1.06 * ones(num_steps, 1) 7.74 * ones(num_steps, 1)];
# define and solve the forward problem
problem = ODEProblem( BoxModel ,c_0 ,tspan ,sources);
sol = solve(problem, alg_hints = "stiff",saveat=time);