Connecting CellML models with ModelingToolkit equations #101
Comments
|
Apologies for the off-topic comment, but perhaps this might be useful in considering this kind of feature in CellMLToolkit.. In work we have been doing with libCellML and code generation with CellML 2.0 models, we generally expect the CellML model to be complete (giving In the very recent work (see cellml/libcellml#1077) the handling of these external variables is moved to the start of the model analysis/code generation. This could begin to support this use-case where such under-defined variables in the CellML model could be detected by the user and set external and so not flagged as an error as it would be taken out of the model analysis. This is an interesting use-case generally that I don't believe we have considered, and I often find myself setting those dummy initial values that I never intend to use and this could solve that problem... |
|
I was able to cook up a working example. using CellMLToolkit, OrdinaryDiffEq, Plots
function custom_stim_model(filename; comp=:membrane, varname=:Istim)
# Load model
doc = CellMLToolkit.load_cellml(filename)
CellMLToolkit.extract_mathml.(doc.xmls)
CellMLToolkit.infer_iv(doc)
# Classify variables as in the original
class = CellMLToolkit.classify_variables(doc)
# We now mark the stimulus as time-dependent
class[CellMLToolkit.Var(comp, varname)] = true
# Construct the systems
systems = CellMLToolkit.subsystems(doc, class)
# Define stimulus
t = CellMLToolkit.get_ivₚ(doc)
stimulus = getproperty(systems[comp], varname) ~ ifelse(t < 1.0, 0.5, 0.0)
# Build the system
odesys = ODESystem(
[connections; stimulus],
t,
systems = collect(values(systems)),
name = CellMLToolkit.gensym(:custom_stim_model)
)
# Simplify
odesys_s = CellMLToolkit.structural_simplify(odesys)
post_sub = CellMLToolkit.post_substitution(doc, systems)
CellMLToolkit.@set! odesys_s.eqs = CellMLToolkit.substitute_eqs(equations(odesys_s), post_sub)
# Update default values
CellMLToolkit.@set! odesys_s.defaults = Dict(CellMLToolkit.find_list_value(doc, vcat(parameters(odesys_s), states(odesys_s))))
return CellModel(doc, odesys_s)
end
br = custom_stim_model("models/beeler_reuter_1977_nostim.cellml.xml")
prob = ODEProblem(br, (0.0, 400.0))
sol = solve(prob, TRBDF2())
plot(sol, vars=[br.sys.membrane₊V])Where the file is simply the provided Beeler-Reuter model without the stimulus protocol component, i.e. models/beeler_reuter_1977_nostim.cellml.xmlThis also works for the forced oscillator above and we can connect multiple models with this approach. However, this example heavily relies on the internals of this package. Hence, I would ask if we could add some of the internals to the public API and keep them at least somewhat stable between versions. I am also not sure if the new equation is constructed as intended ( |
We love PRs! So we can definitely add stuff to the public API This has been weighing on me for a while now: we may want to think about redoing this package to use libcellml with a jll. Then we can merge efforts with the libcellml people instead of having a bunch of slightly different implementations. I tried to do a binary builder thing JuliaPackaging/Yggdrasil#4947 but there was funky stuff in the build that I couldn't quite get right. I know this isn't exactly relevant to the problem you're solving but since you're actively doing CellMLTk stuff, I figured I'd let you know in case you want to pick it up. |
|
Sorry I could not really find time yet to tackle this. Just wanting to check if someone has dumped time into the rework. If not, what is the general idea here? Related to this, my vision would be that we could try to extract components from the CellML file into MTK components. I have especially the idea in mind that we should have the possibility to e.g. extract the individual ion channels in their respective formulation (hodgkin-huxley, markov chain,...) or the calcium model or the sarcomere model, if we can either discover them automatically or if they are baked into the CellML model. |
|
No one has done a substantial rework yet. Extracting this into the component based structure of MTK would be the right thing to do. |
|
In JuliaCon, we talked about identifying CellML components and using Conductor.jl as a bridge to MTK. Of course, the problem is that CellML files do not include the relevant metadata. It should be possible to use pattern-matching to do the job (at least partially); however, it is not trivial and only worth it if there is enough interest. |
|
Indeed, if we have no meta information at all, then pattern matching (and normal forms) become something useful, but I need to think more about the exact patterns here. However, maybe a first step could be to extract components which are actually defined, as e.g. in this cellml tutorial https://models.cellml.org/e/e1/tutorial . Still, not straight forward and I have to check libcellml again on how to do it. |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
First, thanks for the great ecosystem! ModelingToolkit really looks promising.
I think I run into a limitation with the current design of CellMLToolkit. From my understanding it simply parses the cellml document and directly parses and simplifies everything, followed by constructing the ODESystem. This fails if the cellml description is incomplete. However, I think it is a valid use case to connect the missing variables manually. To give a very simple example, let us consider a forced van der Pol oscillator.
Parsing this model fails, because f is not connected to anything:
However, it would be nice to have access to a partially parsed model, such that one could connect the missing components via ModelingToolkit.jl either directly by connecting the partial model with other partial models from other CellML files. To give a rough sketch on what I have in mind:
The text was updated successfully, but these errors were encountered: