Search in sources :

Example 1 with RbmKineticLaw

use of cbit.vcell.model.RbmKineticLaw in project vcell by virtualcell.

the class ReactionRuleKineticsPropertiesPanel method updateKineticChoice.

private void updateKineticChoice() {
    RateLawType newRateLawType = (RateLawType) getKineticsTypeComboBox().getSelectedItem();
    // if same as current kinetics, don't create new one
    if (reactionRule == null || reactionRule.getKineticLaw().getRateLawType().equals(newRateLawType)) {
        return;
    }
    RbmKineticLaw newRbmKineticLaw = new RbmKineticLaw(reactionRule, newRateLawType);
    reactionRule.setKineticLaw(newRbmKineticLaw);
    if (newRateLawType != RateLawType.MassAction) {
        reactionRule.setReversible(false);
        isReversibleCheckBox.setSelected(reactionRule.isReversible());
        isReversibleCheckBox.setEnabled(false);
    } else {
        isReversibleCheckBox.setSelected(reactionRule.isReversible());
        isReversibleCheckBox.setEnabled(true);
    }
}
Also used : RbmKineticLaw(cbit.vcell.model.RbmKineticLaw) RateLawType(cbit.vcell.model.RbmKineticLaw.RateLawType)

Example 2 with RbmKineticLaw

use of cbit.vcell.model.RbmKineticLaw in project vcell by virtualcell.

the class RulebasedTransformer method parseBngOutput.

private void parseBngOutput(SimulationContext simContext, Set<ReactionRule> fromReactions, BNGOutput bngOutput) {
    Model model = simContext.getModel();
    Document bngNFSimXMLDocument = bngOutput.getNFSimXMLDocument();
    bngRootElement = bngNFSimXMLDocument.getRootElement();
    Element modelElement = bngRootElement.getChild("model", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
    Element listOfReactionRulesElement = modelElement.getChild("ListOfReactionRules", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
    List<Element> reactionRuleChildren = new ArrayList<Element>();
    reactionRuleChildren = listOfReactionRulesElement.getChildren("ReactionRule", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
    for (Element reactionRuleElement : reactionRuleChildren) {
        ReactionRuleAnalysisReport rar = new ReactionRuleAnalysisReport();
        boolean bForward = true;
        String rule_id_str = reactionRuleElement.getAttributeValue("id");
        ruleElementMap.put(rule_id_str, reactionRuleElement);
        String rule_name_str = reactionRuleElement.getAttributeValue("name");
        String symmetry_factor_str = reactionRuleElement.getAttributeValue("symmetry_factor");
        Double symmetry_factor_double = Double.parseDouble(symmetry_factor_str);
        ReactionRule rr = null;
        if (rule_name_str.startsWith("_reverse_")) {
            bForward = false;
            rule_name_str = rule_name_str.substring("_reverse_".length());
            rr = model.getRbmModelContainer().getReactionRule(rule_name_str);
            RbmKineticLaw rkl = rr.getKineticLaw();
            try {
                if (symmetry_factor_double != 1.0 && fromReactions.contains(rr)) {
                    Expression expression = rkl.getLocalParameterValue(RbmKineticLawParameterType.MassActionReverseRate);
                    expression = Expression.div(expression, new Expression(symmetry_factor_double));
                    rkl.setLocalParameterValue(RbmKineticLawParameterType.MassActionReverseRate, expression);
                }
            } catch (ExpressionException | PropertyVetoException exc) {
                exc.printStackTrace();
                throw new RuntimeException("Unexpected transform exception: " + exc.getMessage());
            }
            rulesReverseMap.put(rr, rar);
        } else {
            rr = model.getRbmModelContainer().getReactionRule(rule_name_str);
            RbmKineticLaw rkl = rr.getKineticLaw();
            try {
                if (symmetry_factor_double != 1.0 && fromReactions.contains(rr)) {
                    Expression expression = rkl.getLocalParameterValue(RbmKineticLawParameterType.MassActionForwardRate);
                    expression = Expression.div(expression, new Expression(symmetry_factor_double));
                    rkl.setLocalParameterValue(RbmKineticLawParameterType.MassActionForwardRate, expression);
                }
            } catch (ExpressionException | PropertyVetoException exc) {
                exc.printStackTrace();
                throw new RuntimeException("Unexpected transform exception: " + exc.getMessage());
            }
            rulesForwardMap.put(rr, rar);
        }
        rar.symmetryFactor = symmetry_factor_double;
        System.out.println("rule id=" + rule_id_str + ", name=" + rule_name_str + ", symmetry factor=" + symmetry_factor_str);
        keyMap.put(rule_id_str, rr);
        Element listOfReactantPatternsElement = reactionRuleElement.getChild("ListOfReactantPatterns", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
        List<Element> reactantPatternChildren = new ArrayList<Element>();
        reactantPatternChildren = listOfReactantPatternsElement.getChildren("ReactantPattern", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
        for (int i = 0; i < reactantPatternChildren.size(); i++) {
            Element reactantPatternElement = reactantPatternChildren.get(i);
            String pattern_id_str = reactantPatternElement.getAttributeValue("id");
            ReactionRuleParticipant p = bForward ? rr.getReactantPattern(i) : rr.getProductPattern(i);
            SpeciesPattern sp = p.getSpeciesPattern();
            // System.out.println("  reactant id=" + pattern_id_str + ", name=" + sp.toString());
            keyMap.put(pattern_id_str, sp);
            // list of molecules
            extractMolecules(sp, model, reactantPatternElement);
            // list of bonds (not implemented)
            extractBonds(reactantPatternElement);
        }
        Element listOfProductPatternsElement = reactionRuleElement.getChild("ListOfProductPatterns", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
        List<Element> productPatternChildren = new ArrayList<Element>();
        productPatternChildren = listOfProductPatternsElement.getChildren("ProductPattern", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
        for (int i = 0; i < productPatternChildren.size(); i++) {
            Element productPatternElement = productPatternChildren.get(i);
            String pattern_id_str = productPatternElement.getAttributeValue("id");
            ReactionRuleParticipant p = bForward ? rr.getProductPattern(i) : rr.getReactantPattern(i);
            SpeciesPattern sp = p.getSpeciesPattern();
            // System.out.println("  product  id=" + pattern_id_str + ", name=" + sp.toString());
            keyMap.put(pattern_id_str, sp);
            extractMolecules(sp, model, productPatternElement);
            extractBonds(productPatternElement);
        }
        // extract the Map for this rule
        Element listOfMapElement = reactionRuleElement.getChild("Map", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
        List<Element> mapChildren = new ArrayList<Element>();
        mapChildren = listOfMapElement.getChildren("MapItem", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
        for (Element mapElement : mapChildren) {
            String target_id_str = mapElement.getAttributeValue("targetID");
            String source_id_str = mapElement.getAttributeValue("sourceID");
            // System.out.println("Map: target=" + target_id_str + " source=" + source_id_str);
            RbmObject target_object = keyMap.get(target_id_str);
            RbmObject source_object = keyMap.get(source_id_str);
            if (source_object != null) {
                // target_object may be null
                System.out.println("      target=" + target_object + " source=" + source_object);
                Pair<RbmObject, RbmObject> mapEntry = new Pair<RbmObject, RbmObject>(target_object, source_object);
                rar.objmappingList.add(mapEntry);
            }
            Pair<String, String> idmapEntry = new Pair<String, String>(target_id_str, source_id_str);
            rar.idmappingList.add(idmapEntry);
        }
        // ListOfOperations
        Element listOfOperationsElement = reactionRuleElement.getChild("ListOfOperations", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
        List<Element> operationsChildren = new ArrayList<Element>();
        System.out.println("ListOfOperations");
        operationsChildren = listOfOperationsElement.getChildren("StateChange", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
        for (Element operationsElement : operationsChildren) {
            String finalState_str = operationsElement.getAttributeValue("finalState");
            String site_str = operationsElement.getAttributeValue("site");
            RbmObject site_object = keyMap.get(site_str);
            // if(site_object == null) System.out.println("!!! Missing map object " + site_str);
            if (site_object != null) {
                // System.out.println("   finalState=" + finalState_str + " site=" + site_object);
                StateChangeOperation sco = new StateChangeOperation(finalState_str, site_str, site_object);
                rar.operationsList.add(sco);
            }
        }
        System.out.println("AddBond, DeleteBond");
        operationsChildren = listOfOperationsElement.getChildren("AddBond", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
        for (Element operationsElement : operationsChildren) {
            String site1_str = operationsElement.getAttributeValue("site1");
            String site2_str = operationsElement.getAttributeValue("site2");
            RbmObject site1_object = keyMap.get(site1_str);
            RbmObject site2_object = keyMap.get(site2_str);
            // if(site2_object == null) System.out.println("!!! Missing map object " + site2_str);
            if (site1_object != null && site2_object != null) {
                // System.out.println("   site1=" + site1_object + " site2=" + site2_object);
                AddBondOperation abo = new AddBondOperation(site1_str, site2_str, site1_object, site2_object);
                rar.operationsList.add(abo);
            }
        }
        operationsChildren = listOfOperationsElement.getChildren("DeleteBond", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
        for (Element operationsElement : operationsChildren) {
            String site1_str = operationsElement.getAttributeValue("site1");
            String site2_str = operationsElement.getAttributeValue("site2");
            RbmObject site1_object = keyMap.get(site1_str);
            RbmObject site2_object = keyMap.get(site2_str);
            // if(site2_object == null) System.out.println("!!! Missing map object " + site2_str);
            if (site1_object != null && site2_object != null) {
                // System.out.println("   site1=" + site1_object + " site2=" + site2_object);
                DeleteBondOperation dbo = new DeleteBondOperation(site1_str, site2_str, site1_object, site2_object);
                rar.operationsList.add(dbo);
            }
        }
        System.out.println("AddOperation");
        operationsChildren = listOfOperationsElement.getChildren("Add", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
        for (Element operationsElement : operationsChildren) {
            String id_str = operationsElement.getAttributeValue("id");
            RbmObject id_object = keyMap.get(id_str);
            // if(id_object == null) System.out.println("!!! Missing map object " + id_str);
            if (id_object != null) {
                // System.out.println("   id=" + id_str);
                AddOperation ao = new AddOperation(id_str, id_object);
                rar.operationsList.add(ao);
            }
        }
        System.out.println("DeleteMolecules");
        operationsChildren = listOfOperationsElement.getChildren("Delete", Namespace.getNamespace("http://www.sbml.org/sbml/level3"));
        for (Element operationsElement : operationsChildren) {
            String id_str = operationsElement.getAttributeValue("id");
            String delete_molecules_str = operationsElement.getAttributeValue("DeleteMolecules");
            RbmObject id_object = keyMap.get(id_str);
            // if(id_object == null) System.out.println("!!! Missing map object " + id_str);
            int delete_molecules_int = 0;
            if (delete_molecules_str != null) {
                delete_molecules_int = Integer.parseInt(delete_molecules_str);
            }
            if (id_object != null) {
                // System.out.println("   id=" + id_str + ", DeleteMolecules=" + delete_molecules_str);
                DeleteOperation dop = new DeleteOperation(id_str, id_object, delete_molecules_int);
                rar.operationsList.add(dop);
            }
        }
    }
    System.out.println("done parsing xml file");
}
Also used : Element(org.jdom.Element) ArrayList(java.util.ArrayList) Document(org.jdom.Document) ExpressionException(cbit.vcell.parser.ExpressionException) SpeciesPattern(org.vcell.model.rbm.SpeciesPattern) Pair(org.vcell.util.Pair) ReactionRule(cbit.vcell.model.ReactionRule) RbmKineticLaw(cbit.vcell.model.RbmKineticLaw) PropertyVetoException(java.beans.PropertyVetoException) Expression(cbit.vcell.parser.Expression) ReactionRuleParticipant(cbit.vcell.model.ReactionRuleParticipant) RbmObject(org.vcell.model.rbm.RbmObject) Model(cbit.vcell.model.Model)

Example 3 with RbmKineticLaw

use of cbit.vcell.model.RbmKineticLaw in project vcell by virtualcell.

the class RulebasedTransformer method transform.

private void transform(SimulationContext originalSimContext, SimulationContext transformedSimulationContext, ArrayList<ModelEntityMapping> entityMappings, MathMappingCallback mathMappingCallback) throws PropertyVetoException {
    Model newModel = transformedSimulationContext.getModel();
    Model originalModel = originalSimContext.getModel();
    ModelEntityMapping em = null;
    // list of rules created from the reactions; we apply the symmetry factor computed by bionetgen only to these
    Set<ReactionRule> fromReactions = new HashSet<>();
    for (SpeciesContext newSpeciesContext : newModel.getSpeciesContexts()) {
        final SpeciesContext originalSpeciesContext = originalModel.getSpeciesContext(newSpeciesContext.getName());
        // map new and old species contexts
        em = new ModelEntityMapping(originalSpeciesContext, newSpeciesContext);
        entityMappings.add(em);
        if (newSpeciesContext.hasSpeciesPattern()) {
            // it's perfect already and can't be improved
            continue;
        }
        try {
            MolecularType newmt = newModel.getRbmModelContainer().createMolecularType();
            newModel.getRbmModelContainer().addMolecularType(newmt, false);
            MolecularTypePattern newmtp_sc = new MolecularTypePattern(newmt);
            SpeciesPattern newsp_sc = new SpeciesPattern();
            newsp_sc.addMolecularTypePattern(newmtp_sc);
            newSpeciesContext.setSpeciesPattern(newsp_sc);
            RbmObservable newo = new RbmObservable(newModel, "O0_" + newmt.getName() + "_tot", newSpeciesContext.getStructure(), RbmObservable.ObservableType.Molecules);
            MolecularTypePattern newmtp_ob = new MolecularTypePattern(newmt);
            SpeciesPattern newsp_ob = new SpeciesPattern();
            newsp_ob.addMolecularTypePattern(newmtp_ob);
            newo.addSpeciesPattern(newsp_ob);
            newModel.getRbmModelContainer().addObservable(newo);
            // map new observable to old species context
            em = new ModelEntityMapping(originalSpeciesContext, newo);
            entityMappings.add(em);
        } catch (ModelException e) {
            e.printStackTrace();
            throw new RuntimeException("unable to transform species context: " + e.getMessage());
        } catch (PropertyVetoException e) {
            // TODO Auto-generated catch block
            e.printStackTrace();
        }
    }
    ReactionSpec[] reactionSpecs = transformedSimulationContext.getReactionContext().getReactionSpecs();
    for (ReactionSpec reactionSpec : reactionSpecs) {
        if (reactionSpec.isExcluded()) {
            // we create rules only from those reactions which are not excluded
            continue;
        }
        ReactionStep rs = reactionSpec.getReactionStep();
        String name = rs.getName();
        String mangled = TokenMangler.fixTokenStrict(name);
        mangled = newModel.getReactionName(mangled);
        Kinetics k = rs.getKinetics();
        if (!(k instanceof MassActionKinetics)) {
            throw new RuntimeException("Only Mass Action Kinetics supported at this time, reaction \"" + rs.getName() + "\" uses kinetic law type \"" + rs.getKinetics().getName() + "\"");
        }
        boolean bReversible = rs.isReversible();
        ReactionRule rr = new ReactionRule(newModel, mangled, rs.getStructure(), bReversible);
        fromReactions.add(rr);
        MassActionKinetics massActionKinetics = (MassActionKinetics) k;
        List<Reactant> rList = rs.getReactants();
        List<Product> pList = rs.getProducts();
        // counting the stoichiometry - 2A+B means 3 reactants
        int numReactants = 0;
        for (Reactant r : rList) {
            numReactants += r.getStoichiometry();
            if (numReactants > 2) {
                String message = "NFSim doesn't support more than 2 reactants within a reaction: " + name;
                throw new RuntimeException(message);
            }
        }
        int numProducts = 0;
        for (Product p : pList) {
            numProducts += p.getStoichiometry();
            if (bReversible && numProducts > 2) {
                String message = "NFSim doesn't support more than 2 products within a reversible reaction: " + name;
                throw new RuntimeException(message);
            }
        }
        RateLawType rateLawType = RateLawType.MassAction;
        RbmKineticLaw kineticLaw = new RbmKineticLaw(rr, rateLawType);
        try {
            String forwardRateName = massActionKinetics.getForwardRateParameter().getName();
            Expression forwardRateExp = massActionKinetics.getForwardRateParameter().getExpression();
            String reverseRateName = massActionKinetics.getReverseRateParameter().getName();
            Expression reverseRateExp = massActionKinetics.getReverseRateParameter().getExpression();
            LocalParameter fR = kineticLaw.getLocalParameter(RbmKineticLawParameterType.MassActionForwardRate);
            fR.setName(forwardRateName);
            LocalParameter rR = kineticLaw.getLocalParameter(RbmKineticLawParameterType.MassActionReverseRate);
            rR.setName(reverseRateName);
            if (rs.hasReactant()) {
                kineticLaw.setParameterValue(fR, forwardRateExp, true);
            }
            if (rs.hasProduct()) {
                kineticLaw.setParameterValue(rR, reverseRateExp, true);
            }
            // 
            for (KineticsParameter reaction_p : massActionKinetics.getKineticsParameters()) {
                if (reaction_p.getRole() == Kinetics.ROLE_UserDefined) {
                    LocalParameter rule_p = kineticLaw.getLocalParameter(reaction_p.getName());
                    if (rule_p == null) {
                        // 
                        // after lazy parameter creation we didn't find a user-defined rule parameter with this same name.
                        // 
                        // there must be a global symbol with the same name, that the local reaction parameter has overridden.
                        // 
                        ParameterContext.LocalProxyParameter rule_proxy_parameter = null;
                        for (ProxyParameter proxyParameter : kineticLaw.getProxyParameters()) {
                            if (proxyParameter.getName().equals(reaction_p.getName())) {
                                rule_proxy_parameter = (LocalProxyParameter) proxyParameter;
                            }
                        }
                        if (rule_proxy_parameter != null) {
                            // we want to convert to local
                            boolean bConvertToGlobal = false;
                            kineticLaw.convertParameterType(rule_proxy_parameter, bConvertToGlobal);
                        } else {
                            // could find neither local parameter nor proxy parameter
                            throw new RuntimeException("user defined parameter " + reaction_p.getName() + " from reaction " + rs.getName() + " didn't map to a reactionRule parameter");
                        }
                    } else if (rule_p.getRole() == RbmKineticLawParameterType.UserDefined) {
                        kineticLaw.setParameterValue(rule_p, reaction_p.getExpression(), true);
                        rule_p.setUnitDefinition(reaction_p.getUnitDefinition());
                    } else {
                        throw new RuntimeException("user defined parameter " + reaction_p.getName() + " from reaction " + rs.getName() + " mapped to a reactionRule parameter with unexpected role " + rule_p.getRole().getDescription());
                    }
                }
            }
        } catch (ExpressionException e) {
            e.printStackTrace();
            throw new RuntimeException("Problem attempting to set RbmKineticLaw expression: " + e.getMessage());
        }
        rr.setKineticLaw(kineticLaw);
        KineticsParameter[] kpList = k.getKineticsParameters();
        ModelParameter[] mpList = rs.getModel().getModelParameters();
        ModelParameter mp = rs.getModel().getModelParameter(kpList[0].getName());
        ReactionParticipant[] rpList = rs.getReactionParticipants();
        for (ReactionParticipant p : rpList) {
            if (p instanceof Reactant) {
                int stoichiometry = p.getStoichiometry();
                for (int i = 0; i < stoichiometry; i++) {
                    SpeciesPattern speciesPattern = new SpeciesPattern(rs.getModel(), p.getSpeciesContext().getSpeciesPattern());
                    ReactantPattern reactantPattern = new ReactantPattern(speciesPattern, p.getStructure());
                    rr.addReactant(reactantPattern);
                }
            } else if (p instanceof Product) {
                int stoichiometry = p.getStoichiometry();
                for (int i = 0; i < stoichiometry; i++) {
                    SpeciesPattern speciesPattern = new SpeciesPattern(rs.getModel(), p.getSpeciesContext().getSpeciesPattern());
                    ProductPattern productPattern = new ProductPattern(speciesPattern, p.getStructure());
                    rr.addProduct(productPattern);
                }
            }
        }
        // commented code below is probably obsolete, we verify (above) in the reaction the number of participants,
        // no need to do it again in the corresponding rule
        // if(rr.getReactantPatterns().size() > 2) {
        // String message = "NFSim doesn't support more than 2 reactants within a reaction: " + name;
        // throw new RuntimeException(message);
        // }
        // if(rr.getProductPatterns().size() > 2) {
        // String message = "NFSim doesn't support more than 2 products within a reaction: " + name;
        // throw new RuntimeException(message);
        // }
        newModel.removeReactionStep(rs);
        newModel.getRbmModelContainer().addReactionRule(rr);
    }
    for (ReactionRuleSpec rrs : transformedSimulationContext.getReactionContext().getReactionRuleSpecs()) {
        if (rrs == null) {
            continue;
        }
        ReactionRule rr = rrs.getReactionRule();
        if (rrs.isExcluded()) {
            // delete those rules which are disabled (excluded) in the Specifications / Reaction table
            newModel.getRbmModelContainer().removeReactionRule(rr);
            continue;
        }
    }
    // now that we generated the rules we can delete the reaction steps they're coming from
    for (ReactionStep rs : newModel.getReactionSteps()) {
        newModel.removeReactionStep(rs);
    }
    try {
        // we invoke bngl just for the purpose of generating the xml file, which we'll then use to extract the symmetry factor
        generateNetwork(transformedSimulationContext, fromReactions, mathMappingCallback);
    } catch (ClassNotFoundException | IOException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
    }
    System.out.println("Finished RuleBased Transformer.");
}
Also used : Product(cbit.vcell.model.Product) SpeciesContext(cbit.vcell.model.SpeciesContext) Reactant(cbit.vcell.model.Reactant) LocalProxyParameter(cbit.vcell.mapping.ParameterContext.LocalProxyParameter) SpeciesPattern(org.vcell.model.rbm.SpeciesPattern) ExpressionException(cbit.vcell.parser.ExpressionException) KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter) RateLawType(cbit.vcell.model.RbmKineticLaw.RateLawType) HashSet(java.util.HashSet) ReactantPattern(cbit.vcell.model.ReactantPattern) ReactionRule(cbit.vcell.model.ReactionRule) ModelException(cbit.vcell.model.ModelException) ProductPattern(cbit.vcell.model.ProductPattern) RbmObservable(cbit.vcell.model.RbmObservable) RbmKineticLaw(cbit.vcell.model.RbmKineticLaw) IOException(java.io.IOException) MolecularType(org.vcell.model.rbm.MolecularType) PropertyVetoException(java.beans.PropertyVetoException) LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter) ModelParameter(cbit.vcell.model.Model.ModelParameter) ProxyParameter(cbit.vcell.model.ProxyParameter) LocalProxyParameter(cbit.vcell.mapping.ParameterContext.LocalProxyParameter) Expression(cbit.vcell.parser.Expression) ReactionStep(cbit.vcell.model.ReactionStep) Model(cbit.vcell.model.Model) MassActionKinetics(cbit.vcell.model.MassActionKinetics) Kinetics(cbit.vcell.model.Kinetics) MassActionKinetics(cbit.vcell.model.MassActionKinetics) MolecularTypePattern(org.vcell.model.rbm.MolecularTypePattern) ReactionParticipant(cbit.vcell.model.ReactionParticipant)

Example 4 with RbmKineticLaw

use of cbit.vcell.model.RbmKineticLaw in project vcell by virtualcell.

the class RulebasedMathMapping method addParticleJumpProcesses.

private void addParticleJumpProcesses(VariableHash varHash, GeometryClass geometryClass, SubDomain subDomain, HashMap<SpeciesPattern, VolumeParticleSpeciesPattern> speciesPatternMap) throws ExpressionException, MappingException, MathException, PropertyVetoException {
    ArrayList<ReactionRule> rrList = new ArrayList<>();
    for (ReactionRuleSpec rrSpec : getSimulationContext().getReactionContext().getReactionRuleSpecs()) {
        if (!rrSpec.isExcluded()) {
            rrList.add(rrSpec.getReactionRule());
        }
    }
    for (ReactionRule reactionRule : rrList) {
        String jpName = TokenMangler.mangleToSName(reactionRule.getName());
        ArrayList<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
        for (ReactantPattern reactantSpeciesPattern : reactionRule.getReactantPatterns()) {
            reactantParticles.add(speciesPatternMap.get(reactantSpeciesPattern.getSpeciesPattern()));
        }
        ArrayList<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
        for (ProductPattern productSpeciesPattern : reactionRule.getProductPatterns()) {
            productParticles.add(speciesPatternMap.get(productSpeciesPattern.getSpeciesPattern()));
        }
        ArrayList<Action> forwardActions = new ArrayList<Action>();
        ArrayList<Action> reverseActions = new ArrayList<Action>();
        for (ParticleVariable reactant : reactantParticles) {
            forwardActions.add(new Action(reactant, Action.ACTION_DESTROY, new Expression(1.0)));
            reverseActions.add(new Action(reactant, Action.ACTION_CREATE, new Expression(1.0)));
        }
        for (ParticleVariable product : productParticles) {
            forwardActions.add(new Action(product, Action.ACTION_CREATE, new Expression(1.0)));
            reverseActions.add(new Action(product, Action.ACTION_DESTROY, new Expression(1.0)));
        }
        RbmKineticLaw kinetics = reactionRule.getKineticLaw();
        if (kinetics.getRateLawType() == RbmKineticLaw.RateLawType.MassAction) {
            boolean constantMassActionKineticCoefficients = true;
            StringBuffer errorMessage = new StringBuffer();
            Parameter forward_rateParameter = kinetics.getLocalParameter(RbmKineticLawParameterType.MassActionForwardRate);
            Expression substitutedForwardRate = MathUtilities.substituteModelParameters(forward_rateParameter.getExpression(), reactionRule.getNameScope().getScopedSymbolTable());
            if (!substitutedForwardRate.flatten().isNumeric()) {
                errorMessage.append("flattened Kf for reactionRule(" + reactionRule.getName() + ") is not numeric, exp = '" + substitutedForwardRate.flatten().infix() + "'");
                constantMassActionKineticCoefficients = false;
            }
            if (reactionRule.isReversible()) {
                Parameter reverse_rateParameter = kinetics.getLocalParameter(RbmKineticLawParameterType.MassActionReverseRate);
                if (reverse_rateParameter == null || reverse_rateParameter.getExpression() == null) {
                    throw new MappingException("reverse rate constant for reaction rule " + reactionRule.getName() + " is missing");
                }
                Expression substitutedReverseRate = MathUtilities.substituteModelParameters(reverse_rateParameter.getExpression(), reactionRule.getNameScope().getScopedSymbolTable());
                if (!substitutedReverseRate.flatten().isNumeric()) {
                    errorMessage.append("flattened Kr for reactionRule(" + reactionRule.getName() + ") is not numeric, exp = '" + substitutedReverseRate.flatten().infix() + "'");
                    constantMassActionKineticCoefficients = false;
                }
            }
            if (constantMassActionKineticCoefficients) {
                addStrictMassActionParticleJumpProcess(varHash, geometryClass, subDomain, reactionRule, jpName, reactantParticles, productParticles, forwardActions, reverseActions);
            } else {
                throw new MappingException("not mass action: " + errorMessage.toString());
            // addGeneralParticleJumpProcess(varHash, geometryClass, subDomain,
            // reactionRule, jpName,
            // reactantParticles, productParticles,
            // forwardActions, reverseActions);
            }
        } else {
            throw new MappingException("rule-based math generation unsupported for Kinetic Law: " + kinetics.getRateLawType());
        }
    }
// end reactionRules
}
Also used : Action(cbit.vcell.math.Action) ReactionRule(cbit.vcell.model.ReactionRule) ProductPattern(cbit.vcell.model.ProductPattern) VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable) ParticleVariable(cbit.vcell.math.ParticleVariable) ArrayList(java.util.ArrayList) RbmKineticLaw(cbit.vcell.model.RbmKineticLaw) Expression(cbit.vcell.parser.Expression) Parameter(cbit.vcell.model.Parameter) UnresolvedParameter(cbit.vcell.mapping.ParameterContext.UnresolvedParameter) LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter) ModelParameter(cbit.vcell.model.Model.ModelParameter) ReactantPattern(cbit.vcell.model.ReactantPattern)

Example 5 with RbmKineticLaw

use of cbit.vcell.model.RbmKineticLaw in project vcell by virtualcell.

the class XmlReader method getRbmReactionRule.

private ReactionRule getRbmReactionRule(Element reactionRuleElement, Model newModel) throws XmlParseException {
    String n = reactionRuleElement.getAttributeValue(XMLTags.NameAttrTag);
    if (n == null || n.isEmpty()) {
        System.out.println("XMLReader: getRbmReactionRule: name is missing.");
        return null;
    }
    try {
        boolean reversible = Boolean.valueOf(reactionRuleElement.getAttributeValue(XMLTags.RbmReactionRuleReversibleTag));
        // get 1st structure if attribute missing
        String structureName = reactionRuleElement.getAttributeValue(XMLTags.StructureAttrTag, newModel.getStructures()[0].getName());
        Structure structure = newModel.getStructure(structureName);
        ReactionRule reactionRule = new ReactionRule(newModel, n, structure, reversible);
        // we ignore this, name and label are the same thing for now
        String reactionRuleLabel = reactionRuleElement.getAttributeValue(XMLTags.RbmReactionRuleLabelTag);
        // 
        // old style kinetics placed parameter values as attributes
        // look for attributes named ("MassActionKf","MassActionKr","MichaelisMentenKcat","MichaelisMentenKm","SaturableKs","SaturableVmax")
        String[] oldKineticsAttributes = new String[] { XMLTags.RbmMassActionKfAttrTag_DEPRECATED, XMLTags.RbmMassActionKrAttrTag_DEPRECATED, XMLTags.RbmMichaelisMentenKcatAttrTag_DEPRECATED, XMLTags.RbmMichaelisMentenKmAttrTag_DEPRECATED, XMLTags.RbmSaturableKsAttrTag_DEPRECATED, XMLTags.RbmSaturableVmaxAttrTag_DEPRECATED };
        boolean bOldKineticsFound = false;
        for (String oldKineticsAttribute : oldKineticsAttributes) {
            if (reactionRuleElement.getAttribute(oldKineticsAttribute) != null) {
                bOldKineticsFound = true;
            }
        }
        if (bOldKineticsFound) {
            readOldRbmKineticsAttributes(reactionRuleElement, reactionRule);
        } else {
            Element kineticsElement = reactionRuleElement.getChild(XMLTags.KineticsTag, vcNamespace);
            if (kineticsElement != null) {
                String kineticLawTypeString = kineticsElement.getAttributeValue(XMLTags.KineticsTypeAttrTag);
                RbmKineticLaw.RateLawType rateLawType = null;
                if (XMLTags.RbmKineticTypeMassAction.equals(kineticLawTypeString)) {
                    rateLawType = RateLawType.MassAction;
                } else if (XMLTags.RbmKineticTypeMichaelisMenten.equals(kineticLawTypeString)) {
                    rateLawType = RateLawType.MichaelisMenten;
                } else if (XMLTags.RbmKineticTypeSaturable.equals(kineticLawTypeString)) {
                    rateLawType = RateLawType.Saturable;
                } else {
                    throw new RuntimeException("unexpected rate law type " + kineticLawTypeString);
                }
                reactionRule.setKineticLaw(new RbmKineticLaw(reactionRule, rateLawType));
                List<Element> parameterElements = kineticsElement.getChildren(XMLTags.ParameterTag, vcNamespace);
                HashMap<String, ParameterRoleEnum> roleHash = new HashMap<String, ParameterContext.ParameterRoleEnum>();
                roleHash.put(XMLTags.RbmMassActionKfRole, RbmKineticLawParameterType.MassActionForwardRate);
                roleHash.put(XMLTags.RbmMassActionKrRole, RbmKineticLawParameterType.MassActionReverseRate);
                roleHash.put(XMLTags.RbmMichaelisMentenVmaxRole, RbmKineticLawParameterType.MichaelisMentenVmax);
                roleHash.put(XMLTags.RbmMichaelisMentenKmRole, RbmKineticLawParameterType.MichaelisMentenKm);
                roleHash.put(XMLTags.RbmSaturableVmaxRole, RbmKineticLawParameterType.SaturableVmax);
                roleHash.put(XMLTags.RbmSaturableKsRole, RbmKineticLawParameterType.SaturableKs);
                roleHash.put(XMLTags.RbmUserDefinedRole, RbmKineticLawParameterType.UserDefined);
                HashSet<String> xmlRolesToIgnore = new HashSet<String>();
                xmlRolesToIgnore.add(XMLTags.RbmRuleRateRole);
                ParameterContext parameterContext = reactionRule.getKineticLaw().getParameterContext();
                readParameters(parameterElements, parameterContext, roleHash, RbmKineticLawParameterType.UserDefined, xmlRolesToIgnore, newModel);
            }
        }
        Element e1 = reactionRuleElement.getChild(XMLTags.RbmReactantPatternsListTag, vcNamespace);
        getRbmReactantPatternsList(e1, reactionRule, newModel);
        Element e2 = reactionRuleElement.getChild(XMLTags.RbmProductPatternsListTag, vcNamespace);
        getRbmProductPatternsList(e2, reactionRule, newModel);
        reactionRule.checkMatchConsistency();
        return reactionRule;
    } catch (PropertyVetoException | ExpressionException ex) {
        ex.printStackTrace(System.out);
        throw new RuntimeException("failed to parse kinetics for reaction rule '" + n + "': " + ex.getMessage(), ex);
    }
}
Also used : ReactionRule(cbit.vcell.model.ReactionRule) HashMap(java.util.HashMap) Element(org.jdom.Element) RbmKineticLaw(cbit.vcell.model.RbmKineticLaw) ExpressionException(cbit.vcell.parser.ExpressionException) PropertyVetoException(java.beans.PropertyVetoException) ParameterContext(cbit.vcell.mapping.ParameterContext) RateLawType(cbit.vcell.model.RbmKineticLaw.RateLawType) Structure(cbit.vcell.model.Structure) ParameterRoleEnum(cbit.vcell.mapping.ParameterContext.ParameterRoleEnum) HashSet(java.util.HashSet)

Aggregations

RbmKineticLaw (cbit.vcell.model.RbmKineticLaw)9 ReactionRule (cbit.vcell.model.ReactionRule)8 LocalParameter (cbit.vcell.mapping.ParameterContext.LocalParameter)5 Expression (cbit.vcell.parser.Expression)5 ModelParameter (cbit.vcell.model.Model.ModelParameter)4 ProductPattern (cbit.vcell.model.ProductPattern)4 ReactantPattern (cbit.vcell.model.ReactantPattern)4 ArrayList (java.util.ArrayList)4 UnresolvedParameter (cbit.vcell.mapping.ParameterContext.UnresolvedParameter)3 Model (cbit.vcell.model.Model)3 Parameter (cbit.vcell.model.Parameter)3 RateLawType (cbit.vcell.model.RbmKineticLaw.RateLawType)3 SpeciesContext (cbit.vcell.model.SpeciesContext)3 ExpressionException (cbit.vcell.parser.ExpressionException)3 PropertyVetoException (java.beans.PropertyVetoException)3 LocalProxyParameter (cbit.vcell.mapping.ParameterContext.LocalProxyParameter)2 Kinetics (cbit.vcell.model.Kinetics)2 Product (cbit.vcell.model.Product)2 RbmObservable (cbit.vcell.model.RbmObservable)2 Reactant (cbit.vcell.model.Reactant)2