[BUG]: Model extraction cases

[liunelson]

I found some edge cases for template_model_from_sympy_odes.

Case 1

Fixed Birth (l) and proportionate death (m X) processes are interpreted correctly as natural production and natural degradation:

\frac{d S}{d t} = -b S I + l - m S
\frac{d I}{d t} = b S I - g I - m I
\frac{d R}{d t} = g I - m R

However, I had expected the m S terms on d I/d t, d R/d t to be interpreted controlled degradation templates with rate laws m * S.

\frac{d S}{d t} = -b S I + l - m S
\frac{d I}{d t} = b S I - g I - m S
\frac{d R}{d t} = g I - m S

Doing it directly with template_model_from_sympy_odes gives a model that completely ignores all the m S terms:

odes = [
    sympy.Eq(S(t).diff(t), - b * S(t) * I(t) + l - m * S(t)),
    sympy.Eq(I(t).diff(t), b * S(t) * I(t) - g * I(t) - m * S(t)),
    sympy.Eq(R(t).diff(t), g * I(t) - m * S(t)),
]

Doing it within Terarium (which styles the above LaTeX and converts it to SymPy before sending to MIRA) gives a strange model wherein the rate laws are "-I*S*b + l", "I*S*b - I*g", "I*g" (possibly unrelated to MIRA):

• https://app.staging.terarium.ai/projects/a1eb1451-35cf-4f39-814d-d1bd46b08360/workflow/d71c39d9-7bb1-4f23-abf2-fa17d4a19a37?operator=00f41290-365f-42c7-810b-815e73e112b0
• Model (1).json

Case 2

This is a case with "branching ratios":

\frac{d S}{d t} = -b S I
\frac{d I}{d t} = b S I - g I
\frac{d R}{d t} = k g I
\frac{d V}{d t} = (1 - k) g I

This gives a model where I has a natural degradation g I and two controlled productions of R, V, instead of two natural conversions from I into R, V.

If we were to rewrite the 2nd equation as \frac{d I}{d t} = b S I - k g I - (1 - k) g I, then MIRA returns the correct model (where I branches into R, V with ratio k).

Such a branching case was actually involved in a paper from which the UCSD team wanted to extract a model and it required careful reading of the text to realize a rewrite of the equations.

• Equations S1 in the SI of this paper: https://www.sciencedirect.com/science/article/pii/S0264410X21013086#s0120
• Note how I branches into Z, R with ratio eta

Do you have a forthcoming solution to this problem or do you expect users/Terarium to rewrite the equations?

[liunelson]

I have more cases described here:
DARPA-ASKEM/terarium#5805

The issue there (Cases 1, 3) appear to be related to the SymPy's parse_latex adding unnecessary parentheses around multiple terms, causing MIRA to represent the parenthesis'd term to be a single template.

Can simplify_rate_law or MIRA fix this problem or do we need to ditch the SymPy parser and hope that a LLM agent could convert the LaTeX to SymPy better?

Case 3a

odes_latex = [
    r"\frac{d S(t)}{d t} = -b * S(t) * I(t) + l - m * S(t)",
    r"\frac{d I(t)}{d t} = b * S(t) * I(t) - g * I(t) - m * I(t)", 
    r"\frac{d R(t)}{d t} = g * I(t) - m * R(t)",
]

odes_sympy = [
    Eq(Derivative(S(t), t), -m*S(t) + (l + ((-b)*S(t))*I(t))),
    Eq(Derivative(I(t), t), -m*I(t) + (-g*I(t) + (b*S(t))*I(t))),
    Eq(Derivative(R(t), t), g*I(t) - m*R(t))
]

Case 3b

odes_latex = [
    r"\frac{d S(t)}{d t} = -b * S(t) * I(t)",
    r"\frac{d I(t)}{d t} = b * S(t) * I(t) - k * g * I(t) - (1 - k) * g * I(t)", 
    r"\frac{d R(t)}{d t} = k * g * I(t)",
    r"\frac{d V(t)}{d t} = (1 - k) * g * I(t)"
]

odes_sympy = [
    Eq(Derivative(S(t), t), -b*I(t)*S(t)),
    Eq(Derivative(I(t), t), -g*(1 - k)*I(t) + ((b*S(t))*I(t) - g*k*I(t))),
    Eq(Derivative(R(t), t), (g*k)*I(t)),
    Eq(Derivative(V(t), t), (g*(1 - k))*I(t))
]

[bgyori]

For Case 1, I think the term \frac{d R}{d t} = - m S is not physically plausible and doesn't really fit a canonical pattern that could be recognized. A physically plausible controlled degradation would be something like \frac{d R}{d t} = - m R S. Does this come up in practice or is it a hypothetical example?

[bgyori]

For Case 2, I agree this is an ambiguous case (both models produce correct ODEs) and it would be nice to recognize the natural conversions, though not trivial. This requires some further thinking and algorithmic improvement.

[liunelson]

@bgyori
Case 1 was meant to be purely hypothetical - I was added the usual natural death processes (dX/dt = ... - m * X) and tried out this admittedly unphysical variation. Is this feature to ignore unphysical patterns such as this one?

I just updated our LaTeX style guide (which an LLM agent is instructed to follow when cleaning up LaTeX provided by users or upstream service). * will now be used to explicitly denote multiplication to avoid SymPy parse_latex(...) from converting LaTeX a b (1 - g) I(t) to SymPy a * b(t = 1 - g) * I(t).

[liunelson]

I imagine Case 2 is quite nontrivial to tackle automatically, despite how common I see it in model-extraction scenarios. I'm split on whether to try to teach/instruct the equation-styling LLM agent to recognize and expand branching terms. I'll experiment.

[liunelson]

Could you comment on Cases 3a/b? We're trying to figure out how to pass SymPy strings (as opposed to SymPy sympy.core.relational.Equality) to the MIRA function.

Previously, we simply did:

model = template_model_from_sympy_odes([sympy.parsing.latex.parse_latex(ode) for ode in odes_latex])

If MIRA didn't get tripped up by the extra (), we wouldn't have to switch to an LLM solution.

[bgyori]

For Case 3a, I believe we are getting the expected result, despite the parentheses.

ControlledConversion 		 I*S*b
NaturalProduction 		 l
NaturalDegradation 		 S*m
NaturalConversion 		 I*g
NaturalDegradation 		 I*m
NaturalDegradation 		 R*m

In particular, the first two templates look correct in terms of a separate production and conversion template.

[bgyori]

Case 3b appears to be working correctly as well, I get these templates

ControlledConversion 		 I*S*b
NaturalConversion 		 I*g*k
NaturalConversion 		 I*g*(1 - k)

which look correct (just printing some basic details, the actual subjects/objects/controllers are also correct)

[liunelson]

That's quite weird, I get different results from you:

• Case 3a: model.json
  
• Case 3b: model.json
  

I'm also up to the latest MIRA

[liunelson]

Here's the code snippets that I used:

# Case 3a
odes_latex = [
    r"\frac{d S(t)}{d t} = -b * S(t) * I(t) + l - m * S(t)",
    r"\frac{d I(t)}{d t} = b * S(t) * I(t) - g * I(t) - m * I(t)", 
    r"\frac{d R(t)}{d t} = g * I(t) - m * R(t)",
]

odes_sympy = [sympy.parsing.latex.parse_latex(ode) for ode in odes_latex]

__ = [print(ode) for ode in odes_sympy]

model = template_model_from_sympy_odes(odes_sympy)

generate_summary_table(model)

# Case 3b
odes_latex = [
    r"\frac{d S(t)}{d t} = -b * S(t) * I(t)",
    r"\frac{d I(t)}{d t} = b * S(t) * I(t) - k * g * I(t) - (1 - k) * g * I(t)", 
    r"\frac{d R(t)}{d t} = k * g * I(t)",
    r"\frac{d V(t)}{d t} = (1 - k) * g * I(t)"
]

odes_sympy = [sympy.parsing.latex.parse_latex(ode) for ode in odes_latex]

__ = [print(ode) for ode in odes_sympy]

model = template_model_from_sympy_odes(odes_sympy)
generate_summary_table(model)

[liunelson]

@bgyori
Re: SIDARTHE extraction issue

The source of the problem lies with the terms grouped by parentheses.

Here's an example:

[
  "\\frac{d S(t)}{d t} = -S(t) * (\\alpha * I(t) + \\beta * D(t) + \\gamma * A(t) + \\delta * R(t))",
  "\\frac{d I(t)}{d t} = S(t) * (\\alpha * I(t) + \\beta * D(t) + \\gamma * A(t) + \\delta * R(t)) - (\\epsilon + \\zeta + \\lambda) * I(t)",
  "\\frac{d D(t)}{d t} = \\epsilon * I(t) - (\\eta + \\rho) * D(t)",
  "\\frac{d A(t)}{d t} = \\zeta * I(t) - (\\theta + \\mu + \\kappa) * A(t)",
  "\\frac{d R(t)}{d t} = \\eta * D(t) + \\theta * A(t) - (\\nu + \\xi) * R(t)",
  "\\frac{d T(t)}{d t} = \\mu * A(t) + \\nu * R(t) - (\\sigma + \\tau) * T(t)",
  "\\frac{d H(t)}{d t} = \\lambda * I(t) + \\rho * D(t) + \\kappa * A(t) + \\xi * R(t) + \\sigma * T(t)",
  "\\frac{d E(t)}{d t} = \\tau * T(t)"
]

Our agent is instructed to expand all time-dependent algebraic terms enclosed by parentheses (and not for constant terms). However, it is currently only doing so for some terms but not all:

[
  "\\frac{d S(t)}{d t} = -\\alpha * S(t) * I(t) - \\beta * S(t) * D(t) - \\gamma * S(t) * A(t) - \\delta * S(t) * R(t)",
  "\\frac{d I(t)}{d t} = \\alpha * S(t) * I(t) + \\beta * S(t) * D(t) + \\gamma * S(t) * A(t) + \\delta * S(t) * R(t) - (\\epsilon + \\zeta + \\lambda) * I(t)",
  "\\frac{d D(t)}{d t} = \\epsilon * I(t) - (\\eta + \\rho) * D(t)",
  "\\frac{d A(t)}{d t} = \\zeta * I(t) - (\\theta + \\mu + \\kappa) * A(t)",
  "\\frac{d R(t)}{d t} = \\eta * D(t) + \\theta * A(t) - (\\nu + \\xi) * R(t)",
  "\\frac{d T(t)}{d t} = \\mu * A(t) + \\nu * R(t) - (\\sigma + \\tau) * T(t)",
  "\\frac{d H(t)}{d t} = \\lambda * I(t) + \\rho * D(t) + \\kappa * A(t) + \\xi * R(t) + \\sigma * T(t)",
  "\\frac{d E(t)}{d t} = \\tau * T(t)"
]

This generates the equivalent SymPy that retains the parentheses, which seem to trip up MIRA:

import sympy
from sympy import Function, symbols, Eq
 
# Define time variable
t = symbols("t")
 
# Define time-dependent variables
S, I, D, A, R, T, H, E = symbols("S I D A R T H E", cls=Function)
 
# Define constant parameters
alpha, beta, gamma, delta, epsilon, zeta, lambda_, eta, rho, theta, mu, kappa, nu, xi, sigma, tau = symbols("alpha beta gamma delta epsilon zeta lambda eta rho theta mu kappa nu xi sigma tau")
 
# Define the system of differential equations
equation_output = [
    Eq(S(t).diff(t), -alpha * S(t) * I(t) - beta * S(t) * D(t) - gamma * S(t) * A(t) - delta * S(t) * R(t)),
    Eq(I(t).diff(t), alpha * S(t) * I(t) + beta * S(t) * D(t) + gamma * S(t) * A(t) + delta * S(t) * R(t) - (epsilon + zeta + lambda_) * I(t)),
    Eq(D(t).diff(t), epsilon * I(t) - (eta + rho) * D(t)),
    Eq(A(t).diff(t), zeta * I(t) - (theta + mu + kappa) * A(t)),
    Eq(R(t).diff(t), eta * D(t) + theta * A(t) - (nu + xi) * R(t)),
    Eq(T(t).diff(t), mu * A(t) + nu * R(t) - (sigma + tau) * T(t)),
    Eq(H(t).diff(t), lambda_ * I(t) + rho * D(t) + kappa * A(t) + xi * R(t) + sigma * T(t)),
    Eq(E(t).diff(t), tau * T(t))
 ]

This produces a broken model with duplicate templates/transitions (t1, t3, ...):
model.json

Adjusting our agent instructions to impose looser parenthesis expansion gets us the better set of equations:

[
  "\\frac{d S(t)}{d t} = -\\alpha * S(t) * I(t) - \\beta * S(t) * D(t) - \\gamma * S(t) * A(t) - \\delta * S(t) * R(t)",
  "\\frac{d I(t)}{d t} = \\alpha * S(t) * I(t) + \\beta * S(t) * D(t) + \\gamma * S(t) * A(t) + \\delta * S(t) * R(t) - \\epsilon * I(t) - \\zeta * I(t) - \\lambda * I(t)",
  "\\frac{d D(t)}{d t} = \\epsilon * I(t) - \\eta * D(t) - \\rho * D(t)",
  "\\frac{d A(t)}{d t} = \\zeta * I(t) - \\theta * A(t) - \\mu * A(t) - \\kappa * A(t)",
  "\\frac{d R(t)}{d t} = \\eta * D(t) + \\theta * A(t) - \\nu * R(t) - \\xi * R(t)",
  "\\frac{d T(t)}{d t} = \\mu * A(t) + \\nu * R(t) - \\sigma * T(t) - \\tau * T(t)",
  "\\frac{d H(t)}{d t} = \\lambda * I(t) + \\rho * D(t) + \\kappa * A(t) + \\xi * R(t) + \\sigma * T(t)",
  "\\frac{d E(t)}{d t} = \\tau * T(t)"
]

And the correct model:
model_correct.json

However, this instruction change means that it would also expand terms used for scaling (1 - c) and branch ratios (1 - a - b):

a * b * (1 - c) * I(t) -> a * b * I(t) - a * b * c * I(t)
...
= ... - a * I(t) - (1 - a) * I(t) -> - I(t)
= ... + a * I(t)
= ... + (1 - a) * I(t) -> I(t) - a * I(t)

Does MIRA have any particular algo that could tackle this (...) issue in case some groupings slip through?