Hodgkin-Huxley 3: Debugging a model

By the time you have worked through this part of the tutorial you will be able to:

  • Parse an existing CellML file and deserialise it into a model instance;

  • Use the diagnostic Validator class to identify issues in the model’s definition;

  • Use the std::any_cast to retrieve items which need repair from Issue items;

  • Use the Importer class to resolve imports and identify issues; and

  • Use the diagnostic Analyser class to identify issues in the model’s mathematical formulation.

C++ resources

Python resources

CellML resources

Overview

This tutorial looks at how a “broken” model can be investigated and debugged using the diagnostic tools in the three special classes: the Importer, the Analyser, and the Validator.

Step 1: Parse the existing sodium channel model

The Parser class is used to deserialise a CellML string into a Model instance. This means that you’re responsible for finding, opening and reading the *.cellml file into a single string. The parser will then read that string and return a model.

1.a Read a CellML file into a string.

1.b Create a Parser item.

1.c Use the parser to deserialise the contents of the string you’ve read into a model.

1.d Print the parsed model to the terminal for viewing.

Show C++ snippet

    //  1.a 
    //      Read a CellML file into a std::string.
    std::ifstream inFile("sodiumChannelModel_broken.cellml");
    std::stringstream inFileContents;
    inFileContents << inFile.rdbuf();

    //  1.b 
    //      Create a Parser item. 
    auto parser = libcellml::Parser::create();

    //  1.c 
    //      Use the parser to deserialise the contents of the string you've read and return the model.
    auto model = parser->parseModel(inFileContents.str());

    //  1.d 
    //      Print the parsed model to the terminal for viewing.
    printModel(model, false);

Show Python snippet

    #  1.a 
    #      Read a CellML file into a string.
    read_file = open("sodiumChannelModel_broken.cellml")

    #  1.b 
    #      Create a Parser item. 
    parser = Parser()

    #  1.c 
    #      Use the parser to deserialise the contents of the string you've read and return the model.
    model = parser.parseModel(read_file.read())

    #  1.d 
    #      Print the parsed model to the terminal for viewing.
    print_model(model, False)


MODEL: 'SodiumChannelModel'
    UNITS: 5 custom units
        [0]: mV
        [1]: ms
        [2]: per_ms
        [3]: per_mV_ms
        [4]: microA_per_cm2
    COMPONENTS: 4 components
        [0]: controller <--- imported from: 'controller' in 'CircularControllerReference.cellml'
            VARIABLES: 2 variables
                [0]: t
                    └──> sodiumChannel:t [ms]
                [1]: V
                    └──> sodiumChannel:V
        [1]: mGateEquations!
            VARIABLES: 5 variables
                [0]: alpha_m [per_ms]
                [1]: V [mV]
                [2]: beta_m [per_ms]
                [3]: m [dimensionless]
                [4]: t [ms]
        [2]: importedGateM <--- imported from: 'gateEquations' in 'GateModel.cellml'
            VARIABLES: 4 variables
                [0]: alpha_X
                [1]: beta_X
                [2]: X
                [3]: t
        [3]: sodiumChannel
            VARIABLES: 3 variables
                [0]: t [ms]
                    └──> controller:t, sodiumChannelEquations:t [ms]
                [1]: V
                    └──> controller:V, sodiumChannelEquations:V [mV]
                [2]: i_Na [microA_per_cm2]
                    └──> sodiumChannelEquations:i_Na [microA_per_cm2]
            COMPONENT sodiumChannel has 2 child components:
                [0]: sodiumChannelEquations
                    VARIABLES: 8 variables
                        [0]: Na_conductance [mS_per_cm2]
                        [1]: g_Na [mS_per_cm2]
                            └──> sodiumChannelParameters:g_Na [mS_per_cm2]
                        [2]: h [dimensionless]
                            └──> hGate:h [dimensionless]
                        [3]: m [dimensionless]
                            └──> mGate:m [dimensionless]
                        [4]: i_Na [microA_per_cm2]
                            └──> sodiumChannel:i_Na [microA_per_cm2]
                        [5]: V [mV]
                            └──> sodiumChannel:V, hGate:V [mV], mGate:V [mV]
                        [6]: E_Na [mV]
                            └──> sodiumChannelParameters:E_Na [mV]
                        [7]: t [ms]
                            └──> sodiumChannel:t [ms], hGate:t [ms], mGate:t [ms]
                    COMPONENT sodiumChannelEquations has 2 child components:
                        [0]: mGate
                            VARIABLES: 3 variables
                                [0]: t [ms]
                                    └──> sodiumChannelEquations:t [ms]
                                [1]: m [dimensionless]
                                    └──> sodiumChannelEquations:m [dimensionless]
                                [2]: V [mV]
                                    └──> sodiumChannelEquations:V [mV]
                            COMPONENT mGate has 1 child components:
                                [0]: mGateParameters
                                    VARIABLES: 2 variables
                                        [0]: m [dimensionless], initial = 0.05
                                        [1]: i_am_redundant [steradian]
                        [1]: hGate
                            VARIABLES: 3 variables
                                [0]: t [ms]
                                    └──> sodiumChannelEquations:t [ms], hGateEquations:t [ms]
                                [1]: h [dimensionless]
                                    └──> sodiumChannelEquations:h [dimensionless], hGateEquations:h [dimensionless]
                                [2]: V [mV]
                                    └──> sodiumChannelEquations:V [mV], hGateEquations:V [mV]
                            COMPONENT hGate has 2 child components:
                                [0]: hGateParameters
                                    VARIABLES: 1 variables
                                        [0]: h [dimensionless], initial = 0.6
                                [1]: hGateEquations
                                    VARIABLES: 5 variables
                                        [0]: alpha_h [per_ms]
                                            └──> importedGateH:alpha_X
                                        [1]: beta_h [per_ms]
                                            └──> importedGateH:beta_X
                                        [2]: V [mV]
                                            └──> hGate:V [mV]
                                        [3]: h [dimensionless]
                                            └──> hGate:h [dimensionless], importedGateH:X
                                        [4]: t [ms]
                                            └──> hGate:t [ms], importedGateH:t
                                    COMPONENT hGateEquations has 1 child components:
                                        [0]: importedGateH <--- imported from: 'i_dont_exist' in 'GateModel.cellml'
                                            VARIABLES: 4 variables
                                                [0]: alpha_X
                                                    └──> hGateEquations:alpha_h [per_ms]
                                                [1]: beta_X
                                                    └──> hGateEquations:beta_h [per_ms]
                                                [2]: X
                                                    └──> hGateEquations:h [dimensionless]
                                                [3]: t
                                                    └──> hGateEquations:t [ms]
                [1]: sodiumChannelParameters
                    VARIABLES: 2 variables
                        [0]: g_Na [mS_per_cm2], initial = 120
                            └──> sodiumChannelEquations:g_Na [mS_per_cm2]
                        [1]: E_Na [mV]
                            └──> sodiumChannelEquations:E_Na [mV]

Step 2: Validate the parsed model

Create a Validator item and use it to validate the model you’ve just read.

2.a Create a Validator item and validate the model.

Show C++ snippet

    //  2.a 
    //      Create a Validator item and validate the model.
    auto validator = libcellml::Validator::create();
    validator->validateModel(model);

Show Python snippet

    #  2.a 
    #      Create a Validator item and validate the model.
    validator = Validator()
    validator.validateModel(model)

2.b Retrieve any issues from the validator and print them to the terminal.

The validator found 10 issues.
Issue 0: CellML identifiers must not contain any characters other than [a-zA-Z0-9_].
reference: 1.3.1.1
see: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specA03.html?issue=1.3.1.1
stored item type: UNDEFINED

Issue 1: Component 'mGateEquations!' does not have a valid name attribute.
reference: 2.7.1
see: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specB07.html?issue=2.7.1
stored item type: COMPONENT

Issue 2: Variable 'Na_conductance' in component 'sodiumChannelEquations' has a units reference 'mS_per_cm2' which is neither standard nor defined in the parent model.
reference: 2.8.1.2
see: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specB08.html?issue=2.8.1.2
stored item type: VARIABLE

... etc ...

Show C++ snippet

    //  2.b
    //      Retrieve any issues from the validator and print their descriptions and help URL to the terminal.
    //      If you're already familiar with these you can use the printIssues function instead.
    std::cout << "The validator found " << validator->issueCount() << " issues." << std::endl;
    for(size_t i = 0; i < validator->issueCount(); ++i) {
        auto issue = validator->issue(i);
        std::cout << "Issue " << i << ": " << issue->description() << std::endl;
        std::cout << "  reference: "<< issue->referenceHeading() << std::endl;
        std::cout << "  see: " << issue->url() << std::endl;
        std::cout << "  stored item type: " << getCellmlElementTypeFromEnum(issue->cellmlElementType()) << std::endl;
        std::cout << std::endl;
    }

Show Python snippet

    #  2.b
    #      Retrieve any issues from the validator and print their descriptions and help URL to the terminal.
    #      If you're already familiar with these you can use the print_issues function instead.
    print('The validator found {} issues.'.format(validator.issueCount()))
    for i in range(0, validator.issueCount()):
        issue = validator.issue(i)
        print('Issue {}: {}'.format(i,issue.description()))
        print('  reference: {}'.format(issue.referenceHeading()))
        print('  see: {}'.format(issue.url()))
        print('  stored item type: {}'.format(get_cellml_element_type_from_enum(issue.cellmlElementType())))
        print()

Step 3: Repair the parsed model

The messages returned from the validator (and other classes) should (!) have enough information to enable you to know what the problem is. In the case of the validator class, the URL listed contains additional resources and examples related to the issue, as well as its section in the CellML normative specification.

In some situations more than one Issue will be generated from a single cause: this is the case with issues 0 and 1 here:

Issue 0: CellML identifiers must not contain any characters other than [a-zA-Z0-9_].
    reference: 1.3.1.1
    see: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specA03.html?issue=1.3.1.1
Issue 1: Component 'mGateEquations!' does not have a valid name attribute.
    reference: 2.7.1
    see: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specB07.html?issue=2.7.1

3.a Fetch the component using its name, and set the name to something valid. Note that when finding a Component item, setting an optional second argument to true will search the entire encapsulation hierarchy for a component with that name, and not only the direct children of the model. You can follow the URL for information about what makes a valid name, and use the setName function to fix it.

Show C++ snippet

    //  3.a
    //      Fetch the component using its name, and set the name to something valid.  Note that when finding
    //      a Component item, setting an optional second argument to `true` will search the entire encapsulation
    //      hierarchy for a component with that name, and not only the direct children of the model.
    //      You can follow the URL for information about what makes a valid name, and use the Component::setName
    //      function to fix it.
    model->component("mGateEquations!", true)->setName("mGateEquations");

Show Python snippet

    #  3.a
    #      Fetch the component using its name, and set the name to something valid.  Note that when finding
    #      a Component item, setting an optional second argument to `True` will search the entire encapsulation
    #      hierarchy for a component with that name, and not only the direct children of the model.
    #      You can follow the URL for information about what makes a valid name, and use the Component.setName
    #      function to fix it.
    model.component('mGateEquations!', True).setName('mGateEquations')


Issue 2: Variable 'Na_conductance' in component 'sodiumChannelEquations' has a units reference 'mS_per_cm2' which is neither standard nor defined in the parent model.
    reference: 2.8.1.2
    see: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specB08.html?issue=2.8.1.2
Issue 3: Variable 'g_Na' in component 'sodiumChannelEquations' has a units reference 'mS_per_cm2' which is neither standard nor defined in the parent model.
    reference: 2.8.1.2
    see: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specB08.html?issue=2.8.1.2
Issue 4: Variable 'g_Na' in component 'sodiumChannelParameters' has a units reference 'mS_per_cm2' which is neither standard nor defined in the parent model.
    reference: 2.8.1.2
    see: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specB08.html?issue=2.8.1.2

3.b The messages above indicate that we’re missing a Units item named “mS_per_cm2”. Create an appropriate Units item (note that S stands for “siemens”), and add it to your model.

Show C++ snippet

    //  3.b
    //      Create an appropriate Units item (note that S stands for "siemens"), and add it to your model.
    auto mS_per_cm2 = libcellml::Units::create("mS_per_cm2");
    mS_per_cm2->addUnit("siemens", "milli");
    mS_per_cm2->addUnit("metre", "centi", -2);
    model->addUnits(mS_per_cm2);

Show Python snippet

    #  3.b
    #      Create an appropriate Units item (note that S stands for 'siemens'), and add it to your model.
    mS_per_cm2 = Units('mS_per_cm2')
    mS_per_cm2.addUnit('siemens', 'milli')
    mS_per_cm2.addUnit('metre', 'centi', -2)
    model.addUnits(mS_per_cm2)


Issue 5: CellML identifiers must contain one or more basic Latin alphabetic characters.
    reference: 1.3.1.1
    see: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specA03.html?issue=1.3.1.1
Issue 6: Variable 'V' in component 'sodiumChannel' does not have a valid units attribute. The attribute given is ''.
    reference: 2.8.1.2
    see: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specB08.html?issue=2.8.1.2
Issue 8: Variable 'V' in component 'sodiumChannel' has units of '' and an equivalent variable 'V' in component 'sodiumChannelEquations' with non-matching units of 'mV'. The mismatch is:
    reference:
    see:

As with 3.a, here we have more than one issue generated from the same cause: in this case, we haven’t specified units for a variable.

Each issue generated contains a pointer to the item to which it refers. We can retrieve the affected item directly from the issue in one of two ways:

  • retrieving an AnyItem structure (whose “first” attribute is an enum of the CellmlElementType; and “second” attribute is an std::any cast of the item itself); and casting it appropriately, or

  • since we know that the type of item in this error is a VARIABLE, we can call the convenience method Issue::variable() to return the variable which needs attention. (Of course you could retrieve it using the name of its parent component and its name too - this is just another way!)

3.c Check that the item to be returned from the issue is in fact an CellmlElementType::VARIABLE by calling its cellmlElementType() function. Retrieve the variable missing units from the issue. Set its units to be millivolts.

Show C++ snippet

    //  3.c
    //      Check that the item to be returned from the issue is in fact an CellmlElementType::VARIABLE by calling the Issue::type()
    //      function.  Retrieve the variable missing units from the issue.  Set its units to be millivolts.
    auto issue4 = validator->issue(4);
    assert(issue4->cellmlElementType() == libcellml::CellmlElementType::VARIABLE);
    issue4->variable()->setUnits(model->units("mV"));

Show Python snippet

    #  3.c
    #      Check that the item to be returned from the issue is in fact an CellmlElementType.VARIABLE by calling the Issue.type()
    #      function.  Retrieve the variable missing units from the issue.  Set its units to be millivolts.
    issue6 = validator.issue(6)
    print('Issue 6 is a {}'.format(get_cellml_element_type_from_enum(issue6.cellmlElementType())))
    issue6.variable().setUnits(model.units('mV'))

The error below indicates that a child Unit references something which can’t be found.

Issue 7: Units reference 'i_dont_exist' in units 'mV' is not a valid reference to a local units or a standard unit type.
    reference: 2.6.1
    see: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specB06.html?issue=2.6.1

You have a few different options for how to fix this one.

  1. The manual way. The issue description tells us the name of the Units where the issue is, and the name of the units required by one of its children. We can locate both items using these names, remove the problem child, and replace it with a correct one.

  2. The pointer way. As with the example in 3.c, we can retrieve a structure representing the Unit child directly from the issue. Since (as above) we know it will be a Unit item, we can call the issue->unit() function to retrieve it. The Unit item consists of a pointer to its parent Units item, and the index of the relevant child.

  3. The roundabout option. Since the error is saying that units named “i_dont_exist” are missing, we could simply provide them by creating a Units item and adding it to the model.

Units class

3.d Choose your preferred method and use it to retrieve the problem unit attributes and print them all to the terminal. Then fix the issue.

Show C++ snippet

    //  3.d
    //      Choose your preferred method and use it to retrieve the problem unit attributes and print them all to
    //      to the terminal.  Then fix the issue.

    std::string prefix;
    std::string id;
    double exponent;
    double multiplier;
    auto mV = model->units("mV");
    mV->unitAttributes("i_dont_exist", prefix, exponent, multiplier, id);
    std::cout << "The units 'mV' child has attributes: base units = 'i_dont_exist', prefix = '"<< prefix << "', exponent = "<<exponent<<", and multiplier = "<<multiplier <<std::endl;

    // Method 1:
    // mV->removeUnit("i_dont_exist");
    // mV->addUnit("volt", "milli");

    // Method 2:
    auto issue5 = validator->issue(5);
    assert(issue5->cellmlElementType() == libcellml::CellmlElementType::UNIT);
    auto issue5item = issue5->unit();
    issue5item->units()->removeUnit(issue5item->index());
    issue5item->units()->addUnit("volt", "milli");

    // Method 3:
    // auto missingUnits = libcellml::Units::create("i_dont_exist");
    // missingUnits->addUnit("volt", "milli");
    // model->addUnits(missingUnits);

Show Python snippet

    #  3.d
    #      Choose your preferred method and use it to retrieve the problem unit attributes and print them all to
    #      to the terminal.  Then fix the issue.
    prefix = ''
    id = ''
    exponent = 0.0
    multiplier = 0.0
    mV = model.units('mV')
    # prefix, exponent, multiplier, id = mV.unitAttributes('i_dont_exist')
    print(mV.unitAttributes('i_dont_exist'))
    print('The units \'mV\' child has attributes: base units = \'i_dont_exist\', prefix = \'{p}\', exponent = \'{e}\', and multiplier = \'{m}\'.'.format(p=prefix, e=exponent, m=multiplier))

    # Method 1:
    # mV.removeUnit('i_dont_exist')
    # mV.addUnit('volt', 'milli')

    # Method 2:
    issue7 = validator.issue(7)
    print('Issue 7 is a {}'.format(get_cellml_element_type_from_enum(issue7.cellmlElementType())))
    issue7_units = issue7.unit().units()
    issue7_units.removeUnit(issue7.unit().index())
    issue7_units.addUnit('volt', 'milli')

    # Method 3:
    # missing_units = Units('i_dont_exist')
    # missing_units.addUnit('volt', 'milli')
    # model.addUnits(missing_units)

The final validator issue refers to the fact that we need to explicitly specify how other components can access each of the variables in this component.

Issue 9: Variable 't' in component 'sodiumChannelEquations' has no interface type set. The interface type required is 'public_and_private'.
    reference: 3.10.8
    see: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specC10.html?issue=3.10.8

3.e Retrieve the variable either using the issue pointer method, or using the name method, and set its interface to be the required type.

Show C++ snippet

    //  3.e
    //      Retrieve the variable either using the issue pointer method, or using the name method, and set its 
    //      interface to be the required type.
    auto issue7 = validator->issue(7);
    assert(issue7->cellmlElementType() == libcellml::CellmlElementType::VARIABLE);
    issue7->variable()->setInterfaceType("public_and_private");

Show Python snippet

    #  3.e
    #      Retrieve the variable either using the issue pointer method, or using the name method, and set its 
    #      interface to be the required type.
    issue9 = validator.issue(9)
    issue9.variable().setInterfaceType('public_and_private')

3.f Revalidate the model and confirm that the errors have gone.

Model class

  • component (“componentName”, true) will search for the component’s name in the whole of the encapsulation hierarchy.

Component class

Tutorial functions

  • printEncapsulation will output just the names of the components, nested in their encapsulation hierarchy.

3.g Even though the model is free from validation errors, we still need to make sure it represents what we want it to. Print the model to the terminal and check its structure.

3.h Use the addComponent functions to rearrange the components as needed until you have the required structure. Validate the model again.

Show C++ snippet

    //  3.g
    //      Even though the model is free from validation errors, we still need to make sure it
    //      represents what we want it to.  Print the model to the terminal and check its structure.
    printModel(model);

    //  3.h 
    //      Use the addComponent functions to rearrange the components as needed until you have the
    //      required structure.  Validate the model again.             
    auto importedGateM = model->component("importedGateM", true);
    auto mGateEquations = model->component("mGateEquations", true);
    auto mGate = model->component("mGate", true);
    mGateEquations->addComponent(importedGateM);
    mGate->addComponent(mGateEquations);

    validator->validateModel(model);
    printIssues(validator);
    printEncapsulation(model);

Show Python snippet

    #  3.g
    #      Even though the model is free from validation errors, we still need to make sure it
    #      represents what we want it to.  Print the model to the terminal and check its structure.
    print_model(model)

    #  3.h 
    #      Use the addComponent functions to rearrange the components as needed until you have the
    #      required structure.  Validate the model again.             
    imported_gate_m = model.component('importedGateM', True)
    m_gate_equations = model.component('mGateEquations', True)
    m_gate = model.component('mGate', True)
    m_gate_equations.addComponent(imported_gate_m)
    m_gate.addComponent(m_gate_equations)

    validator.validateModel(model)
    print_issues(validator)
    print_encapsulation(model)


Model 'SodiumChannelModel' has 2 components
- Component 'controller' has 0 child components
- Component 'sodiumChannel' has 2 child components
    - Component 'sodiumChannelEquations' has 2 child components
        - Component 'mGate' has 2 child components
            - Component 'mGateParameters' has 0 child components
            - Component 'mGateEquations' has 1 child components
                - Component 'importedGateM' has 0 child components
        - Component 'hGate' has 2 child components
            - Component 'hGateParameters' has 0 child components
            - Component 'hGateEquations' has 1 child components
                - Component 'importedGateH' has 0 child components
    - Component 'sodiumChannelParameters' has 0 child components

Step 4: Resolve the model’s imports

It’s important to remember that the imports are merely instructions for how components or units items should be located: only their syntax is checked by the validator, not that the files exist or contain the required information. To debug the imported aspects of the model, we need to use an Importer class.

To resolve the imports, we need a path to a base location against which any relative file addresses can be resolved. For this tutorial, the files are in the same directory as the code, so simply using an empty string is sufficient.

If they’re another directory, make sure to end your path with a slash, “/”. If they’re in your working directory, enter an empty string.

Importer class

4.a Create an Importer instance and use it to resolve the model.

4.b Similarly to the validator, the importer will log any issues it encounters. Retrieve these and print to the terminal (you can do this manually or using the convenience function as before).

Show C++ snippet

    //  4.a 
    //      Create an Importer instance and use it to resolve the model.
    auto importer = libcellml::Importer::create();
    importer->resolveImports(model, "");

    //  4.b 
    //      Similarly to the validator, the importer will log any issues it encounters.
    //      Retrieve these and print to the terminal (you can do this manually or using the
    //      convenience function as before).
    printIssues(importer);

Show Python snippet

    #  4.a 
    #      Create an Importer instance and use it to resolve the model.
    importer = Importer()
    importer.resolveImports(model, '')

    #  4.b 
    #      Similarly to the validator, the importer will log any issues it encounters.
    #      Retrieve these and print to the terminal (you can do this manually or using the
    #      convenience function as before).
    print_issues(importer)


Recorded 2 issues:

Issue [0] is an ERROR:
    description: Import of component 'importedGateH' from 'GateModel.cellml' requires component named 'i_dont_exist' which cannot be found.
    stored item type: COMPONENT

Fix the issues reported by the importer. This needs to be an iterative process because as more files become available to the importer, the content of those files needs to be checked too.

Component class

4.c We need to change the import reference for the component to be “gateEquations” instead of “i_dont_exist”. You can either retrieve the component using its name or directly from the item stored with the issue.

Show C++ snippet

    //  4.c 
    //      Fix the issues reported by the importer.  This needs to be an iterative process as
    //      more files become available to the importer.
    auto issue0 = importer->issue(0);
    issue0->component()->setImportReference("gateEquations");

Show Python snippet

    #  4.c 
    #      Fix the issues reported by the importer.  This needs to be an iterative process as
    #      more files become available to the importer.
    issue0 = importer.issue(0)
    issue0.component().setImportReference('gateEquations')


Issue [1] is a WARNING:
    description: Cyclic dependencies were found when attempting to resolve components in model 'CircularReferences'. The dependency loop is:
        - component 'importedGateH' is imported from 'i_dont_exist' in 'GateModel.cellml';
        - component 'importedGateM' is imported from 'gateEquations' in 'GateModel.cellml';
        - component 'controller' is imported from 'controller' in 'CircularControllerReference.cellml';
        - component 'controller' is imported from 'controller2' in 'CircularControllerReference2.cellml';
        - component 'controller2' is imported from 'controller' in 'CircularControllerReference.cellml'; and
        - component 'controller' is imported from 'controller2' in 'CircularControllerReference2.cellml'.
    stored item type: UNDEFINED

To fix this, we have two options:

  • to open and repair the file which is actually broken, or

  • to switch the import source in this current model to one which doesn’t have circular imports.

It’s included here to highlight the fact that the Importer class opens and instantiates all required dependencies, and that some of those dependencies may have problems of their own … even issues in files that haven’t (yet) been seen at all by you, the user.

4.d In this example we can change the import of the controller component to have url of ‘SodiumChannelController.cellml’.

4.e Resolve the imports again and check that there are no further issues.

Show C++ snippet

    //  4.d
    //      In this example we can change the import of the controller component to have url of
    //      'SodiumChannelController.cellml'.
    model->component("controller", true)->importSource()->setUrl("SodiumChannelController.cellml");

    //  4.e 
    //      Clear the current issues from the importer using the removeAllIssues function.
    //      Resolve the imports again and check that there are no further issues.
    importer->removeAllIssues();
    importer->resolveImports(model, "");
    printIssues(importer);

Show Python snippet

    #  4.d
    #      In this example we can change the import of the controller component to have url of
    #      'SodiumChannelController.cellml'.
    model.component('controller', True).importSource().setUrl('SodiumChannelController.cellml')

    #  4.e 
    #      Clear all issues from the importer using the removeAllIssues() function.
    #      Resolve the imports again and check that there are no further issues.
    importer.removeAllIssues()
    importer.resolveImports(model, '')
    print_issues(importer)

Step 5: Validate the imported dependencies

At this stage we’ve validated the local model, and we’ve used the Importer class to retrieve all of its import dependencies. These dependencies are stored in the importer’s library, and have not yet been validated or analysed.

Importer class

  • libraryCount returns the number of stored models;

  • library returns the model at the given index or given key string;

  • key returns a key string at the given index;

5.a Use a simple loop to validate each of the models stored in the importer’s library.

Show C++ snippet

    //  5.a 
    //      Use a simple loop to validate each of the models stored in the importer's library.
    for(size_t i = 0; i < importer->libraryCount(); ++i) {
        std::cout << "Imported model at key: " << importer->key(i) << std::endl;
        validator->validateModel(importer->library(i));
        printIssues(validator);
    }

Show Python snippet

    #  5.a 
    #      Use a simple loop to validate each of the models stored in the importer's library.
    for i in range(0, importer.libraryCount()):
        print('Imported model at key: {}'.format(importer.key(i)))
        validator.validateModel(importer.library(i))
        print_issues(validator)
    

Imported model at key: CircularControllerReference.cellml
Recorded 0 issues!

Imported model at key: CircularControllerReference2.cellml
Recorded 0 issues!

Imported model at key: GateModel.cellml
Recorded 0 issues!

Imported model at key: SodiumChannelController.cellml
Recorded 0 issues!

Note that the two files creating the circular import in 4.a are still in the library. To limit ourselves to only those models which are still relevant as the import dependencies of our repaired model, we can iterate through our model’s ImportSource items instead. As soon as the model’s imports have been resolved, all these will point to instantiated models within the importer.

Model class

ImportSource class

5.b Loop through the model’s import source items and print their URLs to the terminal. You’ll notice that these have been used as the keys in the importer library. Check that the importer library’s models are the same as that attached to the import source item.

Show C++ snippet

    //  5.b
    //      Loop through the model's import source items and print their urls to the terminal.
    //      You'll notice that these have been used as the keys in the importer library.
    //      Check that the importer library's models are the same as that attached to the
    //      import source item.
    for(size_t i = 0; i < model->importSourceCount(); ++i) {
        std::cout << "Import source [" << i << "]:" << std::endl;
        std::cout << "     url = " << model->importSource(i)->url() << std::endl;
        std::cout << "     model = " << model->importSource(i)->model()->name() << std::endl;
        std::cout << "     library[url] = " << importer->library(model->importSource(i)->url())->name() << std::endl;
    }

Show Python snippet

    #  5.b
    #      Loop through the model's import source items and print their urls to the terminal.
    #      You'll notice that these have been used as the keys in the importer library.
    #      Check that the importer library's models are the same as that attached to the
    #      import source item.
    for i in range(0, model.importSourceCount()):
        print('Import source [{}]:'.format(i))
        print('     url = {}'.format(model.importSource(i).url()))
        print('     model = {}'.format(model.importSource(i).model()))
        print('     library[url] = {}'.format(importer.library(model.importSource(i).url())))


Import source [0]:
    url = GateModel.cellml
    model = 0x7ff61265b3f0
    library[url] = 0x7ff61265b3f0
Import source [1]:
    url = SodiumChannelController.cellml
    model = 0x7ff6141003c0
    library[url] = 0x7ff6141003c0

Step 6: Analyse the model(s)

As with the validator, the Analyser class is a diagnostic class which will check whether the mathematical representation is ready for simulation. This involves making sure that variables are contained in equations, that integrated variables have initial conditions, and that there are no over- or under-constrained sets of equations. Since this model uses imports, the real mathematical model is hidden from the Analyser (just as it was from the Validator). The way around this is to use the Importer class to create a flat (ie: import-free) version of the same model. If the flat model meets the analyser’s checks, then the importing version will too.

Analyser class

6.a Create an Analyser instance and pass in the model for analysis.

6.b Retrieve and print the issues from the analysis to the screen. We expect to see messages related to un-computed variables, since anything which is imported is missing from this model.

Recorded 19 issues:
Issue [0] is an ERROR:
    description: Variable 'V' in component 'controller' is not computed.
    stored item type: VARIABLE
Issue [1] is an ERROR:
    description: Variable 't' in component 'controller' is not computed.
    stored item type: VARIABLE
Issue [2] is an ERROR:
    description: Variable 'alpha_h' in component 'hGateEquations' is not computed.
    stored item type: VARIABLE

... etc ...

6.c Create a flattened version of the model print it to the screen. Notice that any comments indicating that a component was an import have been removed as these components have been instantiated in the flattened model.

6.d Analyse the flattened model and print the issues to the screen.

Show C++ snippet

    //  6.a 
    //      Create an Analyser instance and pass in the model for analysis.
    //      Useful functions: Analyser::analyseModel(Model)
    auto analyser = libcellml::Analyser::create();
    analyser->analyseModel(model);

    //  6.b 
    //      Retrieve and print the issues from the analysis to the screen.  We expect to see 
    //      messages related to un-computed variables, since anything which is imported is 
    //      missing from this model.
    printIssues(analyser);

    //  6.c 
    //      Create a flattened version of the model print it to the screen.
    //      Notice that any comments indicating that a component was an import have been
    //      removed as these components have been instantiated in the flattened model.
    //      Useful functions:
    //          - Importer::flattenModel(Model) will return a flattened copy.
    auto flatModel = importer->flattenModel(model);
    printModel(flatModel);

    //  6.d 
    //      Analyse the flattened model and print the issues to the screen.
    analyser->analyseModel(flatModel);
    printIssues(analyser);

Show Python snippet

    #  6.a 
    #      Create an Analyser instance and pass in the model for analysis.
    #      Useful functions: Analyser.analyseModel(Model)
    analyser = Analyser()
    analyser.analyseModel(model)

    #  6.b 
    #      Retrieve and print the issues from the analysis to the screen.  We expect to see 
    #      messages related to un-computed variables, since anything which is imported is 
    #      missing from this model.
    print_issues(analyser)

    print()
    #  6.c 
    #      Create a flattened version of the model print it to the screen.
    #      Notice that any comments indicating that a component was an import have been
    #      removed as these components have been instantiated in the flattened model.
    #      Useful functions:
    #          - Importer.flattenModel(Model) will return a flattened copy.
    flat_model = importer.flattenModel(model)
    print_model(flat_model)

    #  6.d 
    #      Analyse the flattened model and print the issues to the screen.
    analyser.analyseModel(flat_model)
    print_issues(analyser)


Issue [0] is an ERROR:
    description: Variable 't' in component 'importedGateM' and variable 't' in component
             'importedGateH' cannot both be the variable of integration.
    stored item type: VARIABLE

The issue returned from the analyser says that we’re trying to use two different variables as the base variable of integration, and the CellML code generation facility (which the analyser is tied to) does not support this yet. It’s still valid CellML though! In this example, the real problem is that these two variables are talking about the same thing, but haven’t been connected to one another yet.

Variable class

6.e Create any necessary variable equivalences so that these two variables are connected. You can refer to your printout of the model’s structure to help if need be, and remember that only variables in a sibling or parent/child relationship can be connected.

6.f Re-flatten and re-analyse the model and print the issues to the terminal.

Show C++ snippet

    //  6.e 
    //      Create any necessary variable equivalences so that these two variables are connected. You
    //      can refer to your printout of the model's structure to help if need be, and remember that only
    //      variables in a sibling or parent/child relationship can be connected.
    //      Useful function: Variable::addEquivalence(v1, v2) will create an equivalence between the 
    //                       variables v1 and v2.

    libcellml::Variable::addEquivalence(model->component("importedGateM", true)->variable("t"), 
                                        model->component("mGateEquations", true)->variable("t"));
    libcellml::Variable::addEquivalence(model->component("mGate", true)->variable("t"), 
                                        model->component("mGateEquations", true)->variable("t"));

    //  6.f Reflatten and re-analyse the model and print the issues to the terminal.
    analyser->analyseModel(importer->flattenModel(model));
    printIssues(analyser);

Show Python snippet

    #  6.e 
    #      Create any necessary variable equivalences so that these two variables are connected. You
    #      can refer to your printout of the model's structure to help if need be, and remember that only
    #      variables in a sibling or parent/child relationship can be connected.
    #      Useful function: Variable.addEquivalence(v1, v2) will create an equivalence between the 
    #                       variables v1 and v2.

    Variable.addEquivalence(model.component('importedGateM', True).variable('t'), 
                                        model.component('mGateEquations', True).variable('t'))
    Variable.addEquivalence(model.component('mGate', True).variable('t'), 
                                        model.component('mGateEquations', True).variable('t'))

    #  6.f Re-flatten and re-analyse the model and print the issues to the terminal.
    analyser.analyseModel(importer.flattenModel(model))
    print_issues(analyser)


Recorded 13 issues:
Issue [0] is an ERROR:
    description: Variable 'X' in component 'importedGateM' is used in an ODE, but it is not initialised.
    stored item type: VARIABLE
Issue [1] is an ERROR:
    description: Variable 'alpha_X' in component 'importedGateM' is not computed.
    stored item type: VARIABLE

Now we see the importance of checking iteratively for issues in the analyser class. The nature of this class means that frequently it is unable to continue processing when an issue is encountered. It’s not unusual to fix one issue only to find twenty more! Two of the errors reported deal with non-initialised variables. Looking at the model printout we can see that this is because the integrated variable “X” (in both the imported gates) hasn’t been connected to its local variable “h” or “m” in the appropriate parameters component.

6.g Create all required connections needed to connect these variables. Re-flatten, re-analyse and print the issues to the terminal.

Show C++ snippet

    //  6.g
    //      Create all required connections needed to connect these variables.
    //      Re-flatten, re-analyse and print the issues to the terminal.
    libcellml::Variable::addEquivalence(model->component("importedGateM", true)->variable("X"), 
                                        model->component("mGateEquations", true)->variable("m"));
    libcellml::Variable::addEquivalence(model->component("mGateParameters", true)->variable("m"), 
                                        model->component("mGateEquations", true)->variable("m"));
    libcellml::Variable::addEquivalence(model->component("importedGateH", true)->variable("X"), 
                                        model->component("hGateEquations", true)->variable("h"));
    libcellml::Variable::addEquivalence(model->component("hGateParameters", true)->variable("h"), 
                                        model->component("hGateEquations", true)->variable("h"));                                    
    analyser->analyseModel(importer->flattenModel(model));
    printIssues(analyser);

Show Python snippet

    #  6.g
    #      Create all required connections needed to connect these variables.
    #      Re-flatten, re-analyse and print the issues to the terminal.
    Variable.addEquivalence(model.component('importedGateM', True).variable('X'), 
                                        model.component('mGateEquations', True).variable('m'))
    Variable.addEquivalence(model.component('mGateParameters', True).variable('m'), 
                                        model.component('mGateEquations', True).variable('m'))
    Variable.addEquivalence(model.component('importedGateH', True).variable('X'), 
                                        model.component('hGateEquations', True).variable('h'))
    Variable.addEquivalence(model.component('hGateParameters', True).variable('h'), 
                                        model.component('hGateEquations', True).variable('h'))                                    
    analyser.analyseModel(importer.flattenModel(model))
    print_issues(analyser)

The nice thing about issues in this class is that frequently a few issues refer to the same single problem. The remainder of the issues reported deal with variables that are not computed. This could mean any one of:

  1. the variable is not included in any equations (it’s completely unused);

  2. the variable is included in an equation, but the equation can’t be evaluated (contains some other un-computed variable(s));

  3. the variable is a constant that should have a value assigned; or

  4. the variable hasn’t been connected to the rest of its definition (usually it’s this one!).

Because the “is not computed” errors are cascading by nature, frequently fixing just one will resolve many others.

  • C++: printEquivalentVariableSet with the variable argument

  • Python: print_equivalent_variable_set with the variable argument

Hints for this tutorial:

  • There is at least one of each kind of problem;

  • There’s a convenience function provided (see below) which will print the equivalent variable set for a given variable. You can use the item stored by each issue and this function to check for missing connections.

  • the addEquivalence function returns a boolean indicating success or otherwise. If you check this as you go it will alert you quickly if you’re trying to connect to a variable that’s not found.

6.h From the printout of your model and the issues listed, determine what needs to happen in order to make the model viable, and do it. Check that your final analysis contains no issues.

Show C++ snippet

    //  6.h
    //      From the printout of your model and the issues listed, determine what needs to happen in 
    //      order to make the model viable, and do it.  Check that your final analysis contains no issues.

    // Connect the mGate to its surroundings.
    libcellml::Variable::addEquivalence(model->component("importedGateM", true)->variable("alpha_X"), 
                                        model->component("mGateEquations", true)->variable("alpha_m"));
    libcellml::Variable::addEquivalence(model->component("importedGateM", true)->variable("beta_X"), 
                                        model->component("mGateEquations", true)->variable("beta_m"));
    libcellml::Variable::addEquivalence(model->component("mGate", true)->variable("V"), 
                                        model->component("mGateEquations", true)->variable("V"));
    libcellml::Variable::addEquivalence(model->component("mGate", true)->variable("m"), 
                                        model->component("mGateEquations", true)->variable("m"));

    // E_Na in sodiumChannelParameters needs to be initialised to 40.
    model->component("sodiumChannelParameters", true)->variable("E_Na")->setInitialValue(40);

    // i_am_redundant in mGateParameters is not required.
    model->component("mGateParameters", true)->removeVariable("i_am_redundant");

    analyser->analyseModel(importer->flattenModel(model));
    printIssues(analyser);

Show Python snippet

    #  6.h
    #      From the printout of your model and the issues listed, determine what needs to happen in 
    #      order to make the model viable, and do it.  Check that your final analysis contains no issues.

    # Connect the m_gate to its surroundings.
    Variable.addEquivalence(model.component('importedGateM', True).variable('alpha_X'), 
                                        model.component('mGateEquations', True).variable('alpha_m'))
    Variable.addEquivalence(model.component('importedGateM', True).variable('beta_X'), 
                                        model.component('mGateEquations', True).variable('beta_m'))
    Variable.addEquivalence(model.component('mGate', True).variable('V'), 
                                        model.component('mGateEquations', True).variable('V'))
    Variable.addEquivalence(model.component('mGate', True).variable('m'), 
                                        model.component('mGateEquations', True).variable('m'))

    # E_Na in sodiumChannelParameters needs to be initialised to 40.
    model.component('sodiumChannelParameters', True).variable('E_Na').setInitialValue(40)

    # i_am_redundant in m_gateParameters is not required.
    model.component('mGateParameters', True).removeVariable('i_am_redundant')

    analyser.analyseModel(importer.flattenModel(model))
    print_issues(analyser)

Step 7: Serialise and print the repaired model

7.a Create a Printer instance and use it to print the CellML-formatted version of the repaired model to a string. Remember we’ll still be printing the original version of the model, not the flattened one!

7.b Write the string to a file named “SodiumChannelModel.cellml”; you will use this in Tutorial 4.

Show C++ snippet

    //  7.a 
    //      Create a Printer instance and use it to print the CellML-formatted version of
    //      the repaired model to a string.  Remember we'll still be printing the original
    //      version of the model, not the flattened one!
    auto printer = libcellml::Printer::create();
    auto modelString = printer->printModel(model);

    //  7.b 
    //      Write the string to a file named "SodiumChannelModel.cellml".
    std::ofstream outFile("SodiumChannelModel.cellml");
    outFile << modelString;
    outFile.close();

Show Python snippet

    #  7.a 
    #      Create a Printer instance and use it to print the CellML-formatted version of
    #      the repaired model to a string.  Remember we'll still be printing the original
    #      version of the model, not the flattened one!
    printer = Printer()
    printed_model_string = printer.printModel(model)
    
    #  7.b 
    #      Write the string to a file named 'SodiumChannelModel.cellml'.
    write_file = open('SodiumChannelModel.cellml', 'w')
    write_file.write(printed_model_string)
    write_file.close()