Recommended language/environment for large scale semi-continuous biological models

Question

We have a fairly large (maybe 1000 equations) differential-algebraic equation model written in ACSLX, an obsolete modelling environment similar to Modelica. The model represents the evolution of a biological system through time, with several forcing functions such as atmospheric temperature which are usually provided on a daily or sub-daily timestep. ACSLX translates the model, including a DAE solver, to C, which is then compiled into a DLL. We run this DLL inside a somewhat complex GUI, often in batch mode consisting of ~100 simulation years. We are in Windows.

We need to convert this model into a new language/environment for future-proofing reasons, and are debating the best option.

• it should be free (to avoid the expense and licensing pain associated with things like Matlab).
• it should be well supported and relatively popular
• the model code should be able to be easily read by biological scientists
• it would be advantageous if it were able to be run within a development environment such as RStudio or Spyder
• it should be able to be called from our GUI (e.g. as a DLL)
• it should be fast
• it should facilitate model calibration, sensitivity analysis, etc.

We currently use C++, Fortran, Javascript, R, SmallTalk for various other projects. Other options seem to be Modelica, Python, C, Java, C#. It's not clear whether we need a strong DE solver like ODEPACK or whether we could write our own, simple time stepper, since the forcing function tends to negate the benefit of efficient ODE solvers.

We have a similar model which is written in Fortran 90 and compiles to a DLL which we run/calibrate from RStudio. This works quite well, although Fortran expertise is hard to come by and it's a little hard to debug. Spyder/Python would allow both sides to be done in a single language, including debugging, although the performance hit might be a problem. Although Cython or Dask might help there. Modelica seems like another option, although support seems patchy, and it might be slow to run.

We would really appreciate any suggestions/thoughts/experience from the SO community.

Thank you for your help.

Answer 1

You should consider giving Julia a try. Let me explain what's going on in the design space right now that would be of interest to you. Full disclosure I am the lead developer of JuliaDiffEq.

JuliaDiffEq and DifferentialEquations.jl has a large feature set dedicated to efficiently integrating computationally-difficult differential equations. It has a simple high level interface in a high-level language of growing popularity, so it has a low barrier to entry. It can be used via development tools like Juno or VSCode. The Julia programming language has all of the mutation features you'd expect in a low level language to optimize your derivative functions to be non-allocating and fast in ways you can't do with higher level "vectorized code". It compiles all of your code (along with doing things like interprocedural optimization and inlining) to get you the speed you'd expect from standard (type-stable) Julia functions.

And the key is that the stack goes very deep. You mentioned you want access to ODEPACK, well this will give you access to a lot of ODE solvers and a lot of DAE solvers which includes Netlib and ODEPACK items like DASKR, LSODA, and Sundials. For example, the components like Sundials.jl are direct auto-generated interfaces directly to the underlying C++ libraries, so you can directly translate any IDA C++ code to Julia code and use that. You can use that to directly call the sensitivity analysis and preconditioner features of Sundials, or use them through the higher level interface. Our ecosystem includes high level features like sensitivity analysis and parameter estimation.

That's if you want to write the differential equations directly. The step above that is modeling features. There are modeling libraries being built on top of the JuliaDiffEq components. We have a small chemical reaction network modeling setup being created

But it's likely not what you're looking for. What you want is something like Modelica, and that is being created as Modia.jl. An overview can be found in this video. That is being built on top of Sundials's IDA and does things like index reduction automatically. In fact, it's being developed by the developers of Modelica as the next way forward. It's not released quite yet since there's a few things we are cleaning up (like compiling the KLU sparse linear solvers into our auto-generated Sundials binaries since these are the most efficient choice for this application), but IMO Hilding's "new Modelica" already looks great so I hope it can be released soon.

Also, Julia is a language that can be compiled. This is at a very early stage, but ahead-of-time compilation (AOT) tools like static-julia and PackageCompiler.jl show this in action.

Lastly, it's all free, has an active community along with chatrooms (and one specifically for differential equations). The team is quite large: we have thirty people on the team, and had 10 mentors for 6 Google Summer of Code projects specifically related to differential equations in Julia. So while young, that demonstrates a lot more developer support and activity than many (most?) other open source differential equation ecosystems. And anything that's missing can be written in Julia (and be fast etc.). I think that hits your points, but if not you can always leave a feature request.

Answer 2

We decided to go with C++, because it

• is free
• is fast
• has a freely available and widely supported numerical integration library (boost/odeint)
• is easy to integrate with our larger simulation driver/GUI
• is widely known
• makes standalone simulation easily available to biological scientists (because they can use RStudio as an interactive simulation environment via Rcpp)

Apart from C++, Python seems to be the other option that might meet these criteria.

C# was rejected because there is no suitable numerical integration library and no interactive simulation environment.

I am translating the CSL code to C++ using R scripts to parse the CSL code. This is proving surprisingly effective. Running the model via Rcpp has also proved straightforward. I have made the project available on Github at

https://github.com/woodwards/csl2cpp