function dudt = circadian_deriv ( t, u )

%*****************************************************************************80
%
%% circadian_deriv() evaluates the right hand side of circadian_ode().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 July 2026
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Albert Goldbeter,
%    A model for circadian oscillation in the Drosophilia period protein (PER),
%    Proceedings of the Royal Society of London, B,
%    Volume 261, pages 319-324, 1995.
%
%  Input:
%
%    real t: the value of the independent variable.
%
%    real u(5): the values of the dependent variables at time T.
%
%  Output:
%
%    real dudt(5): the right hand side of the ODE.
%
  [ kappa1, kappa2, kappas, ...
    k1, k2, k3, k4, kd, ki, km, ...
    n, nud, num, nus, ...
    v1, v2, v3, v4, ... 
    t0, u0, tstop ] = circadian_parameters ( );

  M  = u(1);
  P0 = u(2);
  P1 = u(3);
  P2 = u(4);
  PN = u(5);

  dMdt  =   nus * ki^n / ( ki^n + PN^n ) ...
          - num * M    / ( km   + M );

  dP0dt =   kappas * M ...
          - v1 * P0 / ( k1 + P0 ) ...
          + v2 * P1 / ( k2 + P1 );

  dP1dt =   v1 * P0 / ( k1 + P0 ) ...
          - v2 * P1 / ( k2 + P1 ) ...
          - v3 * P1 / ( k3 + P1 ) ...
          + v4 * P2 / ( k4 + P2 );

  dP2dt =   v3 * P1 / ( k3 + P1 ) ...
          - v4 * P2 / ( k4 + P2 ) ...
          - kappa1 * P2 ...
          + kappa2 * PN ...
          - nud * P2 / ( kd + P2 );

  dPNdt = kappa1 * P2 - kappa2 * PN;

  dudt = [ dMdt; dP0dt; dP1dt; dP2dt; dPNdt ];

  return
end

