ParametricFunction from ParametricNDSolveValue changes when evaluated?

Question

I am new to using ParametricNDSolveValue (and ParametricFunctions in general), so I am hoping that I have made a simple error. First I set up a system of differential equations:

rateeqs = {a'[t] == -k1a a[t] - k12 a[t] c70gs + k21 b[t] c60gs, 
           b'[t] == -k1b b[t] - k21 b[t] c60gs + k12 a[t] c70gs, 
           a[0] == a0, b[0] == b0};
c60gs = c70gs = 5;
maxTime = 0.0025;
e60 = 19060;
e70 = 948;
fitFunc[t_] = e60 a[t] +e70 b[t];
params = {k1a, k1b, k12, k21, a0, b0};
initGuesses = {8000, 100, 4500, 2000, 5. 10^-8, 8 10^-7};

Now I find a solution using ParametricNDSolveValue:

solution = ParametricNDSolveValue[
             Append[rateeqs, y[t] == fitFunc[t]], 
             y, {t, 0, maxTime}, params]

I can plot the results at my initial guesses for the parameters this way:

Plot[(solution @@ initGuesses)[t], {t, 0, maxTime}, PlotRange -> Full]

I can then do a sensitivity analysis on one of the parameters following the model in the Documentation Center:

Plot[Evaluate[((solution @@ params)[t] + 
              {0, 1, -1} 1000 D[(solution @@ params)[t], k1a]) /. 
              MapThread[(#1 -> #2) &, {params, initGuesses}]], 
     {t, 0, maxTime}, Filling -> {2 -> {3}}]

However, when I edit that sensitivity analysis line to look at a different parameter, say k1b, and reexecute it, I get errors:

Switching it back to k1a still gives errors. Trying to reexecute the original plot (without the sensitivity analysis) also now gives errors. Reexecuting the assignment of solution then allows me to replot, and I can do the sensitivity analysis on k1b. But then switching back to k1a introduces the errors.

So something in the sensitivity analysis statement is modifying solution in a way that it can no longer be used. What is going on?

Edited to Add: This issue also seems to be happening inside NonlinearModelFit when I apply it to my ParametricFunction, preventing me from getting to a solution.

Further Edit: It seems to be in the derivative somewhere. If I execute:

D[(solution@@params)[t],k1a]/.MapThread[(#1->#2)&,{params,initGuesses}]

I get an InterpolatingFunction. If I then execute the same line with k1b in place of k1a, I get an error (Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions.). If I redefine solution and then try them in the other order, k1b gives an Interpolatingfunction and k1a gives the error.

It's worse than I thought: If I define solution, and then set solution2=solution, evaluating the derivative of either screws up the other.

Narrowing down the problem: The problem seems to be specifying y[t]==fitFunc[t] in the equations list. It looks like the issue goes away if I leave that out and then make the return expression fitFunc[t]. I am still very curious about what exactly is going wrong here, but I seem to have a workaround.

Answer 1

I think there is a bug here, worth reporting to support.

As a workaround, you can include an ODE for y instead of a constraint:

solution = ParametricNDSolveValue[
    Join[rateeqs, {y'[t] == fitFunc'[t], y[0] == fitFunc[0]}],
    y,
    {t,0,maxTime},
    params
];

Then both:

Plot[
    Evaluate @ ReplaceAll[
        (solution @@ params)[t] + {0, 1, -1} 1000 D[(solution @@ params)[t], k1a],
        MapThread[(#1 -> #2) &, {params, initGuesses}]
    ], 
    {t, 0, maxTime},
    Filling -> {2 -> {3}}
]

and:

Plot[
    Evaluate @ ReplaceAll[
        (solution @@ params)[t] + {0, 1, -1} 1000 D[(solution @@ params)[t], k1b],
        MapThread[(#1 -> #2) &, {params, initGuesses}]
    ], 
    {t, 0, maxTime},
    Filling -> {2 -> {3}}
]

produce the expected plots.