add GSL ODE solver to SLiM

[nobrien97]

Hi!

I'm following up on a discussion from slim-discuss on solving ODEs in SLiM (https://groups.google.com/g/slim-discuss/c/9KrBreuMW6A). My use case is for evolving gene regulatory networks in SLiM (see https://doi.org/10.1371/journal.pgen.1011289), but I imagine that an ODE solver would be helpful for a variety of use cases. I'm happy to contribute towards this if you can include the appropriate GSL code.

There is a problem I foresee that might throw a spanner in the works: ODE solvers solve functions defined at compile-time. One option might be to allow for only particular "forms" of ODE and have the Eidos API dynamically construct the ODE in the solver function itself. This would be quite limited though, and only allow for certain types of ODE to be modelled (e.g. you could use this approach to model gene networks by constructing combinations of Hill functions based on an input adjacency matrix and a vector of parameters, but you are making strong assumptions about how each node interacts with each other in that case). There is also the overhead of constructing the ODE every iteration of the solver, but that might not be too much of a problem (apart from being ugly).

I think the more flexible approach is to code-generate based on a user-defined ODE and runtime compile it. There will likely be problems around cross-platform support, but I think it might still be worth the headache considering the limitations of the other approach. Further, this functionality might be useful for #431 and/or #581 in the future. Happy to hear any other ideas though!

Another question is what sort of syntax do we want for defining the ODE? An intuitive choice might be to define an ODE function in Eidos script. However, a very helpful alternative to include would also be SBML parsing (https://sbml.org/documents/what-is-sbml/). SBML is used in systems biology for storing data on biological networks (e.g. reactants, species, concentrations, kinetic equations/ODEs). It would be very convenient to load a SBML document into SLiM and have it construct an ODE on the fly without needing to predefine it in Eidos script. Another library which parses SBML is roadrunner (this was the frontrunner for me to include in my fork if I had to do it that way).

Regarding evolving the network, there are other considerations. In my model/approach, I treat the coefficients/parameters of the ODE as quantitative traits and solve the ODE for each individual based on their parameter values. Then I apply selection to some measurement of the ODE solution (e.g. in the paper above, the total expression/the AUC). But there are many different measurements you could take - e.g. the time it takes to reach an equilibrium, whether it oscillates or not, maximum expression, etc. So I think the Eidos API should return a dataframe of time points and solution values, so that users can do with that data what they wish in Eidos for flexibility.
I also haven't modelled gene duplications or deletions, which would add or remove nodes from the network and adjust the ODEs. If we go for the code-generation route, that sort of behaviour would be possible.

Anyway, I think that we need to think carefully about this - seems like it isn't a trivial add, but I think it would be quite useful. Keen to hear your thoughts.

Cheers,
Nick

[bhaller]

Hi @nobrien97! Thanks for opening this issue.

Lots to ponder here. I'm not sure roadrunner is a viable option; I don't want SLiM to have any external dependencies (that has been a design stake from the beginning, and is important for ease-of-use and portability), and roadrunner appears to be quite a large piece of software so just incorporating its code into SLiM doesn't seem viable.

I think we should look into the GSL option more. It does not appear to be limited to functions defined at compile time; the doc at https://www.gnu.org/software/gsl/doc/html/ode-initval.html# indicates that the GSL defines a type called gsl_odeiv2_system that "defines a general ODE system with arbitrary parameters", and then their solvers work with that data type to solve it. That probably isn't fast as the compile-time/llvm approach, but... how essential is speed here? I don't imagine you're trying to solve massively complex sets of equations here, right? Are you doing a solve once per tick, or once per individual per tick, or what? I think we ought to look into the GSL avenue at least far enough to decide that it is unacceptable for some reason, such as performance; if we decide that then we can look further, but it is by far the simplest option since the GSL is already a part of SLiM (although some source files that are presently excluded will doubtless need to be added).

SBML is interesting but looks very complex. Again I'm wondering what level of complexity we really need to aim for. I imagine the systems biology folks are using SBML to represent massively complex networks with thousands of equations, but again, you're not doing that in SLiM, right? Maybe it would be easiest to develop some kind of basic, minimal SBML capable of parsing what we expect people to actually want to do in SLiM, and just support that? It looks like there are "levels" of SBML, and level 1 doesn't look too crazy perhaps. But I'd suggest that step one ought to be getting the ODE solver itself working; SBML support is a "nice to have" that could be added later, and it seems possible – even likely – that the best solution would turn out to be an SBML-to-Eidos code generator utility that is a completely separate project from SLiM. (Then it could depend on whatever external libraries proved useful, without forcing SLiM itself to have those external dependencies.)

The alternative, as you say, would be to generate C/C++ code, compile that on demand using things like roadrunner, and dynamically link it in. That's... conceivable, but honestly sounds like it would be a never-ending nightmare of problems on different platforms, and I would certainly not be able to maintain it myself. It sounds like biting off more than we can chew. Is that level of complexity really needed? If we went this route, I'd certainly want to make sure that it is cordoned off such that any platform-specific issues cause problems only for users who actually try to use this feature (as opposed to making the build of SLiM itself fail due to problems in this module).

I take it you're comfortable with programming. So I'd suggest that the next step ought to be for you to experiment with the GSL's ODE solving capabilities, in a separate project just written in C and linking against the regular GSL as an external link dependency. See how well it works, whether it is acceptably fast, whether it successfully converges on solutions or just flails around and errors, etc., and report back. If your report is positive, then I think the path forward would be pretty clear. If not, we'll have to mull it over. Does that sound good? I'm way out of my depth on this stuff (I think I got a D in differential equations, lol), so this will need to be driven largely by you and other interested parties, but I'll help where I can. :->

[nobrien97]

Thanks Ben,

I'm not sure roadrunner is a viable option; I don't want SLiM to have any external dependencies (that has been a design stake from the beginning, and is important for ease-of-use and portability), and roadrunner appears to be quite a large piece of software so just incorporating its code into SLiM doesn't seem viable.

Sure, I agree that roadrunner is much too large to include into SLiM proper. I linked it just to show how others had managed to load and evaluate SBML models at runtime. GSL should be able to do the trick for our purposes, I think.

I think we should look into the GSL option more. It does not appear to be limited to functions defined at compile time; the doc at https://www.gnu.org/software/gsl/doc/html/ode-initval.html# indicates that the GSL defines a type called gsl_odeiv2_system that "defines a general ODE system with arbitrary parameters", and then their solvers work with that data type to solve it.

The gsl_odeiv2_system type contains a function pointer to the system of equations, as well as an optional Jacobian matrix (only required for some solution methods), the dimension of the system of equations, and a pointer to the parameter values. In the example, the function is given by:

int
func (double t, const double y[], double f[],
      void *params)
{
  (void)(t); /* avoid unused parameter warning */
  double mu = *(double *)params;
  f[0] = y[1];                                                  // <- dX/dt
  f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);         // <- dY/dt
  return GSL_SUCCESS;
}

and that is fed into an initializer list for the gsl_odeiv2_system, along with the other parameters:

  gsl_odeiv2_system sys = {func, jac, 2, &mu};

My confusion is around that function pointer - the system of equations needs to be defined in C somewhere, right?

That probably isn't fast as the compile-time/llvm approach, but... how essential is speed here? I don't imagine you're trying to solve massively complex sets of equations here, right? Are you doing a solve once per tick, or once per individual per tick, or what?

I'm solving per individual each tick, and so far only a handful of fairly simple equations per system. So far that has been fine. I'm not worried about performance for GSL - there really aren't huge differences between the C ODE solvers that I've noticed (I've used boost::odeint and Ascent previously, they look pretty similar to GSL).

I imagine the systems biology folks are using SBML to represent massively complex networks with thousands of equations, but again, you're not doing that in SLiM, right?

Yeah, those kinds of systems would probably not be feasible (or even useful for SLiM applications), although I'm not sure how common they are... many of the available models are definitely possible to model - gene regulatory networks, metabolic networks, signalling networks etc. The evolution of these kinds of systems would be very interesting to simulate I think.

But I'd suggest that step one ought to be getting the ODE solver itself working; SBML support is a "nice to have" that could be added later

For sure, I was probably getting ahead of myself - I will focus on the solver first before we consider SBML too much.

The alternative, as you say, would be to generate C/C++ code, compile that on demand using things like roadrunner, and dynamically link it in. That's... conceivable, but honestly sounds like it would be a never-ending nightmare of problems on different platforms, and I would certainly not be able to maintain it myself. It sounds like biting off more than we can chew. Is that level of complexity really needed? If we went this route, I'd certainly want to make sure that it is cordoned off such that any platform-specific issues cause problems only for users who actually try to use this feature (as opposed to making the build of SLiM itself fail due to problems in this module).

I think it is unfortunately necessary if we want to have user-defined ODEs in Eidos, but I don't think we need to solve it with roadrunner. We could perhaps start by only supporting a small feature set for ODE generation, and attempt to keep the generated code as platform and compiler-agnostic as possible (I will have to do some reading into this, I imagine that is much easier said than done). All the library has to do is define the function for the GSL system to reference, it should be pretty small and not require any external dependencies. I imagine the platform-specific problems might come from differences in how dlsym/GetProcAddress work. The downside to this approach is that it requires users to have a compiler on their path - but I figure if you are interested in solving ODEs, you probably have a compiler.

If we only allow ODEs to be defined at the start of simulations (e.g. with an initializeODESystem or similar), we could at least control when/where things can go wrong (and maybe cordon off any bugs into only models which use this system). However, this would prevent updating the ODE structure over time, only the parameters could change (but maybe we should focus on the simplest case first).
I agree that it would be very useful to separate all of this functionality out though from the rest of the code (#581 😉)

I take it you're comfortable with programming. So I'd suggest that the next step ought to be for you to experiment with the GSL's ODE solving capabilities, in a separate project just written in C and linking against the regular GSL as an external link dependency. See how well it works, whether it is acceptably fast, whether it successfully converges on solutions or just flails around and errors, etc., and report back. If your report is positive, then I think the path forward would be pretty clear. If not, we'll have to mull it over. Does that sound good? I'm way out of my depth on this stuff (I think I got a D in differential equations, lol), so this will need to be driven largely by you and other interested parties, but I'll help where I can. :->

Sure, I'll write a little example and come back with a benchmark. Thank you for your time and thoughts on this, I really appreciate it.

[bhaller]

Hi @nobrien97,

I'm not sure roadrunner is a viable option; I don't want SLiM to have any external dependencies (that has been a design stake from the beginning, and is important for ease-of-use and portability), and roadrunner appears to be quite a large piece of software so just incorporating its code into SLiM doesn't seem viable.

Sure, I agree that roadrunner is much too large to include into SLiM proper. I linked it just to show how others had managed to load and evaluate SBML models at runtime. GSL should be able to do the trick for our purposes, I think.

Ah, ok, great.

I think we should look into the GSL option more. It does not appear to be limited to functions defined at compile time; the doc at https://www.gnu.org/software/gsl/doc/html/ode-initval.html# indicates that the GSL defines a type called gsl_odeiv2_system that "defines a general ODE system with arbitrary parameters", and then their solvers work with that data type to solve it.

The gsl_odeiv2_system type contains a function pointer to the system of equations, as well as an optional Jacobian matrix (only required for some solution methods), the dimension of the system of equations, and a pointer to the parameter values. In the example, the function is given by... My confusion is around that function pointer - the system of equations needs to be defined in C somewhere, right?

I see, yes; I didn't look closely enough. Interesting. It is possible that we could design this to actually allow the user to define the ODE equations in Eidos (!). We'd have a single C function that acts as "glue" for that, and a parameter that would indicate which Eidos function ought to be called by the "glue" C function. That would probably be pretty slow, but it'd be cool for sure. :-> Other than that possibility, which might not be fit for purpose, I guess we could either (a) provide a set of predefined functions (with user-selectable parameters) that can be used, or (b) figure out a way to compile the code for the functions at runtime, or in a separately compiled plug-in. Option (a) seems like the way to go for a first pass, for simplicity.

That probably isn't fast as the compile-time/llvm approach, but... how essential is speed here? I don't imagine you're trying to solve massively complex sets of equations here, right? Are you doing a solve once per tick, or once per individual per tick, or what?

I'm solving per individual each tick, and so far only a handful of fairly simple equations per system. So far that has been fine. I'm not worried about performance for GSL - there really aren't huge differences between the C ODE solvers that I've noticed (I've used boost::odeint and Ascent previously, they look pretty similar to GSL).

OK. Yes, I imagine either per-tick or per-individual-per-tick would be the common cases, so speed probably wouldn't be a major factor for most users of this facility.

I imagine the systems biology folks are using SBML to represent massively complex networks with thousands of equations, but again, you're not doing that in SLiM, right?

Yeah, those kinds of systems would probably not be feasible (or even useful for SLiM applications), although I'm not sure how common they are... many of the available models are definitely possible to model... The evolution of these kinds of systems would be very interesting to simulate I think.

Definitely! :->

But I'd suggest that step one ought to be getting the ODE solver itself working; SBML support is a "nice to have" that could be added later

For sure, I was probably getting ahead of myself - I will focus on the solver first before we consider SBML too much.

OK, good. Good to have it in the back of our minds, though.

The alternative, as you say, would be to generate C/C++ code, compile that on demand using things like roadrunner, and dynamically link it in...

I think it is unfortunately necessary if we want to have user-defined ODEs in Eidos, but I don't think we need to solve it with roadrunner. We could perhaps start by only supporting a small feature set for ODE generation, and attempt to keep the generated code as platform and compiler-agnostic as possible (I will have to do some reading into this, I imagine that is much easier said than done). All the library has to do is define the function for the GSL system to reference, it should be pretty small and not require any external dependencies. I imagine the platform-specific problems might come from differences in how dlsym/GetProcAddress work. The downside to this approach is that it requires users to have a compiler on their path - but I figure if you are interested in solving ODEs, you probably have a compiler.

Yes. But as I wrote above, I think the first step is to get this working with some set of predefined built-in functions that the user can reference in their Eidos script without needing any dynamic linking at all. We can just provide a linear function, a quadratic function, etc., as built-ins with parameters that the user can specify, right? Some small set of such functions would probably satisfy a large fraction of what people will want to do; and if a particular function is not supported, people can request it and we can add it. Extending this to user-written C code that is dynamically compiled or linked seems like a stretch goal that we can push out until we've got a basic design working. I want to see what this actually looks like, in terms of an Eidos API and usage in an actual SLiM model, before we overcomplicate things. :->

If we only allow ODEs to be defined at the start of simulations (e.g. with an initializeODESystem or similar), we could at least control when/where things can go wrong (and maybe cordon off any bugs into only models which use this system). However, this would prevent updating the ODE structure over time, only the parameters could change (but maybe we should focus on the simplest case first). I agree that it would be very useful to separate all of this functionality out though from the rest of the code (#581 😉)

I think probably we can allow it to be more dynamic than that. I'm envisioning maybe having the user create a Dictionary or DataFrame object that represents the set of equations, and then having a single Eidos function, solveODESystem() or some such name, that you pass that object to. It solves it using the GSL and returns results to you, probably again as a Dictionary or DataFrame of some kind. If the user wants to solve a given system of equations multiple times with different parameter values, they can just keep the Dictionary/DataFrame around in a global and tweak it before each call to solveODESystem(). I'm hoping that we don't end up needing to add a new object class to Eidos for this; it seems like Dictionary should probably suffice.

I take it you're comfortable with programming. So I'd suggest that the next step ought to be for you to experiment with the GSL's ODE solving capabilities, in a separate project just written in C and linking against the regular GSL as an external link dependency. See how well it works, whether it is acceptably fast, whether it successfully converges on solutions or just flails around and errors, etc., and report back. If your report is positive, then I think the path forward would be pretty clear. If not, we'll have to mull it over. Does that sound good? I'm way out of my depth on this stuff (I think I got a D in differential equations, lol), so this will need to be driven largely by you and other interested parties, but I'll help where I can. :->

Sure, I'll write a little example and come back with a benchmark. Thank you for your time and thoughts on this, I really appreciate it.

Sure. It's an interesting idea, and I've been wanting to extend SLiM to consider gene/regulatory networks for a while now but haven't been sure how to proceed! Thanks for thinking about how to do it in a general-purpose way; these kinds of contributions are invaluable, since you know much more about this area than I do! :->

[petrelharp]

@attrna and I have also implemented GRN models kinda like this, but using simple matrix iteration, which allows for adding in noise. This functionality would be pretty attractive, but coming up with an API seems daunting.