What are the sensitivity equations that can be integrated to find these derivatives in the Hessian matrix? In the paper Gutenkunst et al. 2007, part of the mathematical model analysis involves computing the Hessian matrix
$$H_{j,k}^{\chi^2}=\frac{\text d^2\chi^2}{\text d\log\theta_j\,\text d\log\theta_k}$$
where $j$ and $k$ are indices to denote different model parameters $\theta$, and $\chi^2$, which quantifies the change in model behavior as parameters $\theta$ are varied from their published values $\theta^*$, is given by
$$\chi^2(\theta)=\frac{1}{2N_cN_s}\sum_{s,c}\frac{1}{T_c}\int_0^{T_c}\left[\frac{y_{s,c}(\theta,t)-y_{s,c}(\theta^*,t)}{\sigma_s}\right]^2\,\text dt$$
Where $s$ and $c$ (with corresponding numbers $N_s$ and $N_c$) just denote various species / conditions the model considers (irrelevant for this question), $T_c$ is length of the time interval in the data of condition $c$, $\sigma_s$ is a normalization factor, $y_{s,c}(\theta,t)$ is the model output.
The appendix of the paper combines the two above equations to evaluate $H^{\chi^2}$ at $\theta^*$
$$H_{j,k}^{\chi^2}=\frac{1}{N_cN_s}\sum\frac{1}{T_c\sigma_s^2}\int_0^{T_c}\frac{\text dy_{s,c}(\theta^*,t)}{\text d\log\theta_j}\frac{\text dy_{s,c}(\theta^*,t)}{\text d\log\theta_k}\,\text dt$$
The paper then states that this "is convenient because the first derivatives $\text dy_{s,c}(\theta^*,t)/ \text d \log \theta_j$ can be calculated by integrating sensitivity equations. This avoids the use of finite-difference derivatives, which are troublesome in sloppy systems." I have been trying to search around, but I cannot seem to find anything about these sensitivity equations they mention.
What are the sensitivity equations they mention here that can be integrated to determine the first derivatives in computing the Hessian matrix?
 A: As outlined in the paper by Dickinson and Gelinas, we are trying to find the sensitivity for species $s$ to changes in parameter $\theta_j$ given by
$$Z_{s,j}=\frac{\partial y_s}{\partial\theta_j}$$
If we have our system of ODE's given by
$$\dot y_s=f_s(y_1,y_2,\dots,y_n,t,\boldsymbol\theta)$$
We can determine time ODE's for the sensitivities switching ordering of partial differentiation and using the chain rule:
\begin{align}
\dot Z_{s,j}&=\frac{\partial}{\partial t}\left(\frac{\partial y_s}{\partial\theta_j}\right)\\
&\\
&=\frac{\partial}{\partial \theta_j}\left(\frac{\partial y_s}{\partial t}\right)\\
&\\
&=\frac{\partial}{\partial \theta_j}\left(f_s(\mathbf y,t,\boldsymbol\theta)\right)\\
&\\
&=\frac{\partial f_s}{\partial\theta_j}+\sum_{i=1}^n\frac{\partial f_s}{\partial y_i}\frac{\partial y_i}{\partial\theta_j}\\
&\\
&=\frac{\partial f_s}{\partial\theta_j}+\sum_{i=1}^n\frac{\partial f_s}{\partial y_i}Z_{i,j}\\
\end{align}
Therefore, the determination of the sensitivities $Z_{s,j}=\partial y_s/\partial\theta_j$ involves time integration of the model equations along with the system of differential equations
$$\dot{\mathbf{Z}}_{j}=\frac{\partial \mathbf f}{\partial\theta_j}+\mathbf J\mathbf Z_j$$
where $\mathbf Z_j=[Z_{1,j},Z_{2,j},\dots,Z_{n,j}]$ is the vector of sensitivities for each species corresponding to parameter $\theta_j$, $\mathbf f=[f_1,f_2,\dots,f_n]$ is the vector of the right hand sides of the system of model differential equations, and $\mathbf J$ is the $n\times n$ Jacobian matrix $[\mathbf J]_{s,i}=\partial f_s/\partial y_i$.
For initial conditions, if parameter $\theta_j$ is not an initial condition for $y_s$, then $\dot Z_{s,j}(t=0)=0$. If parameter $\theta_j$ is an initial condition for $y_s$, then $\dot Z_{s,j}(t=0)=1$.
