Examples with SimpleReaction cbit.vcell.model.SimpleReaction used on opensource projects

Example 1 with SimpleReaction

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

the class SBMLImporter method addReactions.

/**
 * addReactions:
 */
protected void addReactions(VCMetaData metaData) {
    if (sbmlModel == null) {
        throw new SBMLImportException("SBML model is NULL");
    }
    ListOf<Reaction> reactions = sbmlModel.getListOfReactions();
    final int numReactions = reactions.size();
    if (numReactions == 0) {
        lg.info("No Reactions");
        return;
    }
    // all reactions
    ArrayList<ReactionStep> vcReactionList = new ArrayList<>();
    // just the fast ones
    ArrayList<ReactionStep> fastReactionList = new ArrayList<>();
    Model vcModel = vcBioModel.getSimulationContext(0).getModel();
    ModelUnitSystem vcModelUnitSystem = vcModel.getUnitSystem();
    SpeciesContext[] vcSpeciesContexts = vcModel.getSpeciesContexts();
    try {
        for (Reaction sbmlRxn : reactions) {
            ReactionStep vcReaction = null;
            String rxnName = sbmlRxn.getId();
            boolean bReversible = true;
            if (sbmlRxn.isSetReversible()) {
                bReversible = sbmlRxn.getReversible();
            }
            // Check of reaction annotation is present; if so, does it have
            // an embedded element (flux or simpleRxn).
            // Create a fluxReaction or simpleReaction accordingly.
            Element sbmlImportRelatedElement = sbmlAnnotationUtil.readVCellSpecificAnnotation(sbmlRxn);
            Structure reactionStructure = getReactionStructure(sbmlRxn, vcSpeciesContexts, sbmlImportRelatedElement);
            if (sbmlImportRelatedElement != null) {
                Element embeddedRxnElement = getEmbeddedElementInAnnotation(sbmlImportRelatedElement, REACTION);
                if (embeddedRxnElement != null) {
                    if (embeddedRxnElement.getName().equals(XMLTags.FluxStepTag)) {
                        // If embedded element is a flux reaction, set flux
                        // reaction's strucure, flux carrier, physicsOption
                        // from the element attributes.
                        String structName = embeddedRxnElement.getAttributeValue(XMLTags.StructureAttrTag);
                        CastInfo<Membrane> ci = SBMLHelper.getTypedStructure(Membrane.class, vcModel, structName);
                        if (!ci.isGood()) {
                            throw new SBMLImportException("Appears that the flux reaction is occuring on " + ci.actualName() + ", not a membrane.");
                        }
                        vcReaction = new FluxReaction(vcModel, ci.get(), null, rxnName, bReversible);
                        vcReaction.setModel(vcModel);
                        // Set the fluxOption on the flux reaction based on
                        // whether it is molecular, molecular & electrical,
                        // electrical.
                        String fluxOptionStr = embeddedRxnElement.getAttributeValue(XMLTags.FluxOptionAttrTag);
                        if (fluxOptionStr.equals(XMLTags.FluxOptionMolecularOnly)) {
                            ((FluxReaction) vcReaction).setPhysicsOptions(ReactionStep.PHYSICS_MOLECULAR_ONLY);
                        } else if (fluxOptionStr.equals(XMLTags.FluxOptionMolecularAndElectrical)) {
                            ((FluxReaction) vcReaction).setPhysicsOptions(ReactionStep.PHYSICS_MOLECULAR_AND_ELECTRICAL);
                        } else if (fluxOptionStr.equals(XMLTags.FluxOptionElectricalOnly)) {
                            ((FluxReaction) vcReaction).setPhysicsOptions(ReactionStep.PHYSICS_ELECTRICAL_ONLY);
                        } else {
                            localIssueList.add(new Issue(vcReaction, issueContext, IssueCategory.SBMLImport_Reaction, "Unknown FluxOption : " + fluxOptionStr + " for SBML reaction : " + rxnName, Issue.SEVERITY_WARNING));
                        // logger.sendMessage(VCLogger.Priority.MediumPriority,
                        // VCLogger.ErrorType.ReactionError,
                        // "Unknown FluxOption : " + fluxOptionStr +
                        // " for SBML reaction : " + rxnName);
                        }
                    } else if (embeddedRxnElement.getName().equals(XMLTags.SimpleReactionTag)) {
                        // if embedded element is a simple reaction, set
                        // simple reaction's structure from element
                        // attributes
                        vcReaction = new SimpleReaction(vcModel, reactionStructure, rxnName, bReversible);
                    }
                } else {
                    vcReaction = new SimpleReaction(vcModel, reactionStructure, rxnName, bReversible);
                }
            } else {
                vcReaction = new SimpleReaction(vcModel, reactionStructure, rxnName, bReversible);
            }
            // set annotations and notes on vcReactions[i]
            sbmlAnnotationUtil.readAnnotation(vcReaction, sbmlRxn);
            sbmlAnnotationUtil.readNotes(vcReaction, sbmlRxn);
            // the limit on the reactionName length.
            if (rxnName.length() > 64) {
                String freeTextAnnotation = metaData.getFreeTextAnnotation(vcReaction);
                if (freeTextAnnotation == null) {
                    freeTextAnnotation = "";
                }
                StringBuffer oldRxnAnnotation = new StringBuffer(freeTextAnnotation);
                oldRxnAnnotation.append("\n\n" + rxnName);
                metaData.setFreeTextAnnotation(vcReaction, oldRxnAnnotation.toString());
            }
            // Now add the reactants, products, modifiers as specified by
            // the sbmlRxn
            addReactionParticipants(sbmlRxn, vcReaction);
            KineticLaw kLaw = sbmlRxn.getKineticLaw();
            Kinetics kinetics = null;
            if (kLaw != null) {
                // Convert the formula from kineticLaw into MathML and then
                // to an expression (infix) to be used in VCell kinetics
                ASTNode sbmlRateMath = kLaw.getMath();
                Expression kLawRateExpr = getExpressionFromFormula(sbmlRateMath);
                Expression vcRateExpression = new Expression(kLawRateExpr);
                // modifier (catalyst) to the reaction.
                for (int k = 0; k < vcSpeciesContexts.length; k++) {
                    if (vcRateExpression.hasSymbol(vcSpeciesContexts[k].getName())) {
                        if ((vcReaction.getReactant(vcSpeciesContexts[k].getName()) == null) && (vcReaction.getProduct(vcSpeciesContexts[k].getName()) == null) && (vcReaction.getCatalyst(vcSpeciesContexts[k].getName()) == null)) {
                            // This means that the speciesContext is not a
                            // reactant, product or modifier : it has to be
                            // added to the VC Rxn as a catalyst
                            vcReaction.addCatalyst(vcSpeciesContexts[k]);
                        }
                    }
                }
                // set kinetics on VCell reaction
                if (bSpatial) {
                    // if spatial SBML ('isSpatial' attribute set), create
                    // DistributedKinetics)
                    SpatialReactionPlugin ssrplugin = (SpatialReactionPlugin) sbmlRxn.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
                    // 'spatial'
                    if (ssrplugin != null && ssrplugin.getIsLocal()) {
                        kinetics = new GeneralKinetics(vcReaction);
                    } else {
                        kinetics = new GeneralLumpedKinetics(vcReaction);
                    }
                } else {
                    kinetics = new GeneralLumpedKinetics(vcReaction);
                }
                // set kinetics on vcReaction
                vcReaction.setKinetics(kinetics);
                // If the name of the rate parameter has been changed by
                // user, or matches with global/local param,
                // it has to be changed.
                resolveRxnParameterNameConflicts(sbmlRxn, kinetics, sbmlImportRelatedElement);
                /**
                 * Now, based on the kinetic law expression, see if the rate
                 * is expressed in concentration/time or substance/time : If
                 * the compartment_id of the compartment corresponding to
                 * the structure in which the reaction takes place occurs in
                 * the rate law expression, it is in concentration/time;
                 * divide it by the compartment size and bring in the rate
                 * law as 'Distributed' kinetics. If not, the rate law is in
                 * substance/time; bring it in (as is) as 'Lumped' kinetics.
                 */
                ListOf<LocalParameter> localParameters = kLaw.getListOfLocalParameters();
                for (LocalParameter p : localParameters) {
                    String paramName = p.getId();
                    KineticsParameter kineticsParameter = kinetics.getKineticsParameter(paramName);
                    if (kineticsParameter == null) {
                        // add unresolved for now to prevent errors in kinetics.setParameterValue(kp,vcRateExpression) below
                        kinetics.addUnresolvedParameter(paramName);
                    }
                }
                KineticsParameter kp = kinetics.getAuthoritativeParameter();
                if (lg.isDebugEnabled()) {
                    lg.debug("Setting " + kp.getName() + ":  " + vcRateExpression.infix());
                }
                kinetics.setParameterValue(kp, vcRateExpression);
                // If there are any global parameters used in the kinetics,
                // and if they have species,
                // check if the species are already reactionParticipants in
                // the reaction. If not, add them as catalysts.
                KineticsProxyParameter[] kpps = kinetics.getProxyParameters();
                for (int j = 0; j < kpps.length; j++) {
                    if (kpps[j].getTarget() instanceof ModelParameter) {
                        ModelParameter mp = (ModelParameter) kpps[j].getTarget();
                        HashSet<String> refSpeciesNameHash = new HashSet<String>();
                        getReferencedSpeciesInExpr(mp.getExpression(), refSpeciesNameHash);
                        java.util.Iterator<String> refSpIterator = refSpeciesNameHash.iterator();
                        while (refSpIterator.hasNext()) {
                            String spName = refSpIterator.next();
                            org.sbml.jsbml.Species sp = sbmlModel.getSpecies(spName);
                            ArrayList<ReactionParticipant> rpArray = getVCReactionParticipantsFromSymbol(vcReaction, sp.getId());
                            if (rpArray == null || rpArray.size() == 0) {
                                // This means that the speciesContext is not
                                // a reactant, product or modifier : it has
                                // to be added as a catalyst
                                vcReaction.addCatalyst(vcModel.getSpeciesContext(sp.getId()));
                            }
                        }
                    }
                }
                // model - local params cannot be defined by rules.
                for (LocalParameter param : localParameters) {
                    String paramName = param.getId();
                    Expression exp = new Expression(param.getValue());
                    String unitString = param.getUnits();
                    VCUnitDefinition paramUnit = sbmlUnitIdentifierHash.get(unitString);
                    if (paramUnit == null) {
                        paramUnit = vcModelUnitSystem.getInstance_TBD();
                    }
                    // check if sbml local param is in kinetic params list;
                    // if so, add its value.
                    boolean lpSet = false;
                    KineticsParameter kineticsParameter = kinetics.getKineticsParameter(paramName);
                    if (kineticsParameter != null) {
                        if (lg.isDebugEnabled()) {
                            lg.debug("Setting local " + kineticsParameter.getName() + ":  " + exp.infix());
                        }
                        kineticsParameter.setExpression(exp);
                        kineticsParameter.setUnitDefinition(paramUnit);
                        lpSet = true;
                    } else {
                        UnresolvedParameter ur = kinetics.getUnresolvedParameter(paramName);
                        if (ur != null) {
                            kinetics.addUserDefinedKineticsParameter(paramName, exp, paramUnit);
                            lpSet = true;
                        }
                    }
                    if (!lpSet) {
                        // check if it is a proxy parameter (specifically,
                        // speciesContext or model parameter (structureSize
                        // too)).
                        KineticsProxyParameter kpp = kinetics.getProxyParameter(paramName);
                        // and units to local param values
                        if (kpp != null && kpp.getTarget() instanceof ModelParameter) {
                            kinetics.convertParameterType(kpp, false);
                            kineticsParameter = kinetics.getKineticsParameter(paramName);
                            kinetics.setParameterValue(kineticsParameter, exp);
                            kineticsParameter.setUnitDefinition(paramUnit);
                        }
                    }
                }
            } else {
                // sbmlKLaw was null, so creating a GeneralKinetics with 0.0
                // as rate.
                kinetics = new GeneralKinetics(vcReaction);
            }
            // end - if-else KLaw != null
            // set the reaction kinetics, and add reaction to the vcell
            // model.
            kinetics.resolveUndefinedUnits();
            // System.out.println("ADDED SBML REACTION : \"" + rxnName +
            // "\" to VCModel");
            vcReactionList.add(vcReaction);
            if (sbmlRxn.isSetFast() && sbmlRxn.getFast()) {
                fastReactionList.add(vcReaction);
            }
        }
        // end - for vcReactions
        ReactionStep[] array = vcReactionList.toArray(new ReactionStep[vcReactionList.size()]);
        vcModel.setReactionSteps(array);
        final ReactionContext rc = vcBioModel.getSimulationContext(0).getReactionContext();
        for (ReactionStep frs : fastReactionList) {
            final ReactionSpec rs = rc.getReactionSpec(frs);
            rs.setReactionMapping(ReactionSpec.FAST);
        }
    } catch (ModelPropertyVetoException mpve) {
        throw new SBMLImportException(mpve.getMessage(), mpve);
    } catch (Exception e1) {
        e1.printStackTrace(System.out);
        throw new SBMLImportException(e1.getMessage(), e1);
    }
}

Example 2 with SimpleReaction

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

the class RbmNetworkGenerator method generateModel.

public static void generateModel(BioModel bioModel, String netfile) throws Exception {
    Model model = bioModel.getModel();
    Map<String, SpeciesContext> speciesMap = new HashMap<String, SpeciesContext>();
    Map<String, ReactionStep> reactionMap = new HashMap<String, ReactionStep>();
    List<ReactionLine> reactionLineList = new ArrayList<ReactionLine>();
    BufferedReader br = new BufferedReader(new StringReader(netfile));
    int reversibleCount = 0;
    int reactionCount = 0;
    while (true) {
        String line = br.readLine();
        if (line == null) {
            break;
        }
        line = line.trim();
        if (line.equals(BEGIN_PARAMETERS)) {
            while (true) {
                String line2 = br.readLine();
                line2 = line2.trim();
                if (line2.length() == 0) {
                    continue;
                }
                if (line2.equals(END_PARAMETERS)) {
                    break;
                }
                StringTokenizer st = new StringTokenizer(line2);
                String token1 = st.nextToken();
                String token2 = st.nextToken();
                String token3 = st.nextToken();
                ModelParameter mp = model.new ModelParameter(token2, new Expression(token3), Model.ROLE_UserDefined, bioModel.getModel().getUnitSystem().getInstance_TBD());
                model.addModelParameter(mp);
            }
        } else if (line.equals(BEGIN_SPECIES)) {
            while (true) {
                String line2 = br.readLine();
                line2 = line2.trim();
                if (line2.length() == 0) {
                    continue;
                }
                if (line2.equals(END_SPECIES)) {
                    break;
                }
                StringTokenizer st = new StringTokenizer(line2);
                // no
                String token1 = st.nextToken();
                // pattern
                String token2 = st.nextToken();
                // initial condition
                String token3 = st.nextToken();
                String newname = token2.replaceAll("\\.", "_");
                newname = newname.replaceAll("[\\(,][a-zA-Z]\\w*", "");
                newname = newname.replaceAll("~|!\\d*", "");
                newname = newname.replaceAll("\\(\\)", "");
                newname = newname.replaceAll("\\)", "");
                SpeciesContext sc = model.createSpeciesContext(model.getStructure(0));
                sc.setName(newname);
                bioModel.getVCMetaData().setFreeTextAnnotation(sc, token2);
                bioModel.getVCMetaData().setFreeTextAnnotation(sc.getSpecies(), token2);
                speciesMap.put(token1, sc);
            }
        } else if (line.equals(BEGIN_REACTIONS)) {
            while (true) {
                String line2 = br.readLine();
                line2 = line2.trim();
                if (line2.length() == 0) {
                    continue;
                }
                if (line2.equals(END_REACTIONS)) {
                    break;
                }
                ++reactionCount;
                StringTokenizer st = new StringTokenizer(line2);
                String token1 = st.nextToken();
                // reactants
                String token2 = st.nextToken();
                // products
                String token3 = st.nextToken();
                // rate
                String token4 = st.nextToken();
                String token5 = st.nextToken();
                boolean bFoundReversible = false;
                Expression rate = new Expression(token4);
                for (ReactionLine rl : reactionLineList) {
                    if (token2.equals(rl.products) && token3.equals(rl.reactants) && token5.equals(rl.ruleLabel + "r")) {
                        ReactionStep rs = reactionMap.get(rl.no);
                        ((MassActionKinetics) rs.getKinetics()).getReverseRateParameter().setExpression(rate);
                        reactionLineList.remove(rl);
                        bFoundReversible = true;
                        break;
                    }
                }
                if (bFoundReversible) {
                    ++reversibleCount;
                    continue;
                }
                ReactionLine rl = new ReactionLine(token1, token2, token3, token5);
                reactionLineList.add(rl);
                SimpleReaction reaction = model.createSimpleReaction(model.getStructure(0));
                reactionMap.put(token1, reaction);
                reaction.setModel(model);
                bioModel.getVCMetaData().setFreeTextAnnotation(reaction, line2);
                MassActionKinetics kinetics = new MassActionKinetics(reaction);
                reaction.setKinetics(kinetics);
                st = new StringTokenizer(token2, ",");
                while (st.hasMoreTokens()) {
                    String t = st.nextToken();
                    SpeciesContext sc = speciesMap.get(t);
                    if (sc != null) {
                        boolean bExists = false;
                        for (ReactionParticipant rp : reaction.getReactionParticipants()) {
                            if (rp instanceof Reactant && rp.getSpeciesContext() == sc) {
                                rp.setStoichiometry(rp.getStoichiometry() + 1);
                                bExists = true;
                                break;
                            }
                        }
                        if (!bExists) {
                            reaction.addReactant(sc, 1);
                        }
                    }
                }
                st = new StringTokenizer(token3, ",");
                while (st.hasMoreTokens()) {
                    String t = st.nextToken();
                    SpeciesContext sc = speciesMap.get(t);
                    if (sc != null) {
                        boolean bExists = false;
                        for (ReactionParticipant rp : reaction.getReactionParticipants()) {
                            if (rp instanceof Product && rp.getSpeciesContext() == sc) {
                                rp.setStoichiometry(rp.getStoichiometry() + 1);
                                bExists = true;
                                break;
                            }
                        }
                        if (!bExists) {
                            reaction.addProduct(sc, 1);
                        }
                    }
                }
                kinetics.getForwardRateParameter().setExpression(rate);
            }
        }
    }
    System.out.println(model.getNumSpecies() + " species added");
    System.out.println(model.getNumReactions() + " reactions added");
    System.out.println(reversibleCount + " reversible reactions found");
    if (reactionCount != model.getNumReactions() + reversibleCount) {
        throw new RuntimeException("Reactions are not imported correctly!");
    }
}

Example 3 with SimpleReaction

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

the class SBPAXKineticsExtractor method extractKineticsExactMatch.

public static boolean extractKineticsExactMatch(ReactionStep reaction, Process process) throws ExpressionException, PropertyVetoException {
    // we try a perfect match first based on the existence of a SBOTerm in the kinetic law
    if (process.getControl() == null) {
        // no control means no kinetic law - nothing more to do
        return true;
    }
    Control control = process.getControl();
    ArrayList<SBEntity> sbEntities = control.getSBSubEntity();
    for (SBEntity sbE : sbEntities) {
        // the only SBSubEntities allowed in a control are kinetic laws
        if (sbE.getID().contains("kineticLaw")) {
            // build list of the parameters of this kinetic law in sboParams
            // params of this kinetic law
            ArrayList<SBOParam> sboParams = new ArrayList<SBOParam>();
            ArrayList<SBEntity> subEntities = sbE.getSBSubEntity();
            for (SBEntity subEntity : subEntities) {
                if (subEntity instanceof SBMeasurable) {
                    SBMeasurable m = (SBMeasurable) subEntity;
                    if (!m.hasTerm()) {
                        // we don't know what to do with a measurable with no SBTerm
                        break;
                    }
                    String termName = m.extractSBOTermAsString();
                    SBOTerm sboT = SBOListEx.sboMap.get(termName);
                    System.out.println("    " + sboT.getIndex() + "  " + sboT.getName());
                    SBOParam sboParam = matchSBOParam(sboT);
                    if (m.hasNumber()) {
                        double number = m.getNumber().get(0);
                        sboParam.setNumber(number);
                    }
                    if (m.hasUnit()) {
                        String unit = m.extractSBOUnitAsString();
                        sboParam.setUnit(unit);
                    }
                    sboParams.add(sboParam);
                }
            }
            // find if a kinetic law type exists and if not guesstimate one based on parameters
            // simple rule: if we have a Km param it's MM, otherwise it's mass action
            ArrayList<SBVocabulary> sbTerms = sbE.getSBTerm();
            if (sbTerms.isEmpty()) {
                // return false;
                SBVocabulary sbTerm = new SBVocabulary();
                ArrayList<String> termNames = new ArrayList<String>();
                String id;
                SBOParam kMichaelis = extractMichaelisForwardParam(sboParams);
                if (kMichaelis == null) {
                    // mass action rate law
                    id = new String("SBO:0000012");
                } else {
                    // irreversible Henri-Michaelis-Menten rate law
                    id = new String("SBO:0000029");
                }
                // termNames.add(id);
                // sbTerm.setTerm(termNames);
                sbTerm.setID(id);
                sbTerms.add(sbTerm);
            }
            System.out.println(" kinetic law ID: " + sbTerms.get(0).getID());
            // identify the kinetic law type (mass action, michaelis menten, etc) and bring it in vCell
            for (SBVocabulary sbv : sbTerms) {
                // use for loop even though we only expect 1 SBTerm
                // SBVocabulary id, used to retrieve the kinetic law type
                String vocabularyID = sbv.getID();
                SBOTerm sboT = SBOUtil.getSBOTermFromVocabularyId(vocabularyID);
                System.out.println(vocabularyID + "   " + sboT.getName());
                SBOParam kForward;
                SBOParam kCat;
                SBOParam vM;
                SBOParam kReverse;
                SBOParam kMichaelis;
                Kinetics kinetics;
                MappedKinetics current = matchSBOKineticLaw(sboT);
                switch(current) {
                    case SBO_0000069:
                    case SBO_0000432:
                        // some kinetic laws unknown to vCell will fall through to this category
                        // honestly i don't know what to do with them
                        System.out.println("GeneralKinetics");
                        // TODO: what to do here?
                        return true;
                    case SBO_0000012:
                    case SBO_0000078:
                        System.out.println("MassActionKinetics - reversible");
                        kForward = extractKForwardParam(sboParams);
                        kReverse = extractKReverseParam(sboParams);
                        kinetics = new MassActionKinetics(reaction);
                        reaction.setKinetics(kinetics);
                        setKForwardParam(reaction, kForward, kinetics);
                        setKReverseParam(reaction, kReverse, kinetics);
                        return true;
                    case SBO_0000043:
                        System.out.println("MassActionKinetics - zeroth order irreversible, Kr <- 0  ");
                        kForward = extractKForwardParam(sboParams);
                        kinetics = new MassActionKinetics(reaction);
                        // TODO: what to do here?
                        return true;
                    case SBO_0000044:
                        System.out.println("MassActionKinetics - first order irreversible, Kr <- 0  ");
                        kForward = extractKForwardParam(sboParams);
                        kinetics = new MassActionKinetics(reaction);
                        reaction.setKinetics(kinetics);
                        setKForwardParam(reaction, kForward, kinetics);
                        return true;
                    case SBO_0000028:
                    case SBO_0000029:
                        System.out.println("HMM_IRRKinetics");
                        // TODO: make kCat global variable, set its number and unit in annotation
                        // TODO: make vM global variable, set its number and unit in annotation
                        // get the numbers, if present (may be null)
                        kMichaelis = extractMichaelisForwardParam(sboParams);
                        vM = extractVMForwardParam(sboParams, process);
                        kCat = extractKCatForwardParam(sboParams);
                        kinetics = new HMM_IRRKinetics((SimpleReaction) reaction);
                        try {
                            // TODO: create expression only if kCat != null, otherwise use vM directly (if != null) otherwise ???
                            kinetics.reading(true);
                            setVMForwardParamAsExpression(reaction, kCat, kinetics);
                        } finally {
                            kinetics.reading(false);
                        }
                        reaction.setKinetics(kinetics);
                        setMichaelisForwardParam(reaction, kMichaelis, kinetics);
                        return true;
                    case SBO_0000438:
                        System.out.println("HMMREVKinetics");
                        kinetics = new HMM_REVKinetics((SimpleReaction) reaction);
                        return true;
                    default:
                        // TODO: guessing happens above - if we have nothing by now we need to raise runtime exception
                        // change the code below !!!
                        System.out.println("Unable to match the SBOTerm with any compatible kinetic law.");
                        // found unmapped kinetic law, we'll try to guess a match
                        return false;
                }
            }
        }
    }
    // no SBTerm found so we cannot know for sure the kinetic law, we'll have to guess it
    return false;
}

Example 4 with SimpleReaction

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

the class StochMathMapping method addJumpProcesses.

private void addJumpProcesses(VariableHash varHash, GeometryClass geometryClass, SubDomain subDomain) throws ExpressionException, ModelException, MappingException, MathException {
    // set up jump processes
    // get all the reactions from simulation context
    // ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();---need to take a look here!
    ModelUnitSystem modelUnitSystem = getSimulationContext().getModel().getUnitSystem();
    ReactionSpec[] reactionSpecs = getSimulationContext().getReactionContext().getReactionSpecs();
    for (ReactionSpec reactionSpec : reactionSpecs) {
        if (reactionSpec.isExcluded()) {
            continue;
        }
        // get the reaction
        ReactionStep reactionStep = reactionSpec.getReactionStep();
        Kinetics kinetics = reactionStep.getKinetics();
        // probability parameter from modelUnitSystem
        VCUnitDefinition probabilityParamUnit = modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getTimeUnit());
        // Different ways to deal with simple reactions and flux reactions
        if (// simple reactions
        reactionStep instanceof SimpleReaction) {
            // check the reaction rate law to see if we need to decompose a reaction(reversible) into two jump processes.
            // rate constants are important in calculating the probability rate.
            // for Mass Action, we use KForward and KReverse,
            // for General Kinetics we parse reaction rate J to see if it is in Mass Action form.
            Expression forwardRate = null;
            Expression reverseRate = null;
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction) || kinetics.getKineticsDescription().equals(KineticsDescription.General)) {
                Expression rateExp = new Expression(kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate), reactionStep.getNameScope());
                Parameter forwardRateParameter = null;
                Parameter reverseRateParameter = null;
                if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                    forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward);
                    reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse);
                }
                MassActionSolver.MassActionFunction maFunc = MassActionSolver.solveMassAction(forwardRateParameter, reverseRateParameter, rateExp, reactionStep);
                if (maFunc.getForwardRate() == null && maFunc.getReverseRate() == null) {
                    throw new MappingException("Cannot generate stochastic math mapping for the reaction:" + reactionStep.getName() + "\nLooking for the rate function according to the form of k1*Reactant1^Stoir1*Reactant2^Stoir2...-k2*Product1^Stoip1*Product2^Stoip2.");
                } else {
                    if (maFunc.getForwardRate() != null) {
                        forwardRate = maFunc.getForwardRate();
                    }
                    if (maFunc.getReverseRate() != null) {
                        reverseRate = maFunc.getReverseRate();
                    }
                }
            } else // if it's macro/microscopic kinetics, we'll have them set up as reactions with only forward rate.
            if (kinetics.getKineticsDescription().equals(KineticsDescription.Macroscopic_irreversible) || kinetics.getKineticsDescription().equals(KineticsDescription.Microscopic_irreversible)) {
                Expression Kon = getIdentifierSubstitutions(new Expression(reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_KOn), getNameScope()), reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_Binding_Radius).getUnitDefinition(), geometryClass);
                if (Kon != null) {
                    Expression KonCopy = new Expression(Kon);
                    try {
                        MassActionSolver.substituteParameters(KonCopy, true).evaluateConstant();
                        forwardRate = new Expression(Kon);
                    } catch (ExpressionException e) {
                        throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Problem with Kon parameter in " + reactionStep.getName() + ":  '" + KonCopy.infix() + "', " + e.getMessage()));
                    }
                } else {
                    throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Kon parameter of " + reactionStep.getName() + " is null."));
                }
            }
            boolean isForwardRatePresent = false;
            boolean isReverseRatePresent = false;
            if (forwardRate != null) {
                isForwardRatePresent = true;
            }
            if (reverseRate != null) {
                isReverseRatePresent = true;
            }
            // we process it as forward reaction
            if ((isForwardRatePresent)) /*|| ((forwardRate == null) && (reverseRate == null))*/
            {
                // get jump process name
                String jpName = TokenMangler.mangleToSName(reactionStep.getName());
                // get probability
                Expression exp = null;
                // reactions are of mass action form
                exp = getProbabilityRate(reactionStep, forwardRate, true);
                ProbabilityParameter probParm = null;
                try {
                    probParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, exp, PARAMETER_ROLE_P, probabilityParamUnit, reactionStep);
                } catch (PropertyVetoException pve) {
                    pve.printStackTrace();
                    throw new MappingException(pve.getMessage());
                }
                // add probability to function or constant
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(probParm, geometryClass), getIdentifierSubstitutions(exp, probabilityParamUnit, geometryClass), geometryClass));
                JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, geometryClass)));
                // actions
                ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
                for (int j = 0; j < reacPart.length; j++) {
                    Action action = null;
                    SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[j].getSpeciesContext());
                    if (reacPart[j] instanceof Reactant) {
                        // check if the reactant is a constant. If the species is a constant, there will be no action taken on this species
                        if (// not a constant
                        !simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
                            int stoi = ((Reactant) reacPart[j]).getStoichiometry();
                            action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-stoi));
                            jp.addAction(action);
                        }
                    } else if (reacPart[j] instanceof Product) {
                        // check if the product is a constant. If the product is a constant, there will be no action taken on this species
                        if (// not a constant
                        !simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
                            int stoi = ((Product) reacPart[j]).getStoichiometry();
                            action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(stoi));
                            jp.addAction(action);
                        }
                    }
                }
                // add jump process to compartment subDomain
                subDomain.addJumpProcess(jp);
            }
            if (// one more jump process for a reversible reaction
            isReverseRatePresent) {
                // get jump process name
                String jpName = TokenMangler.mangleToSName(reactionStep.getName()) + PARAMETER_PROBABILITY_RATE_REVERSE_SUFFIX;
                Expression exp = null;
                // reactions are mass actions
                exp = getProbabilityRate(reactionStep, reverseRate, false);
                ProbabilityParameter probRevParm = null;
                try {
                    probRevParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, exp, PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionStep);
                } catch (PropertyVetoException pve) {
                    pve.printStackTrace();
                    throw new MappingException(pve.getMessage());
                }
                // add probability to function or constant
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, geometryClass), getIdentifierSubstitutions(exp, probabilityParamUnit, geometryClass), geometryClass));
                JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, geometryClass)));
                // actions
                ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
                for (int j = 0; j < reacPart.length; j++) {
                    Action action = null;
                    SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[j].getSpeciesContext());
                    if (reacPart[j] instanceof Reactant) {
                        // check if the reactant is a constant. If the species is a constant, there will be no action taken on this species
                        if (// not a constant
                        !simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
                            int stoi = ((Reactant) reacPart[j]).getStoichiometry();
                            action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(stoi));
                            jp.addAction(action);
                        }
                    } else if (reacPart[j] instanceof Product) {
                        // check if the product is a constant. If the product is a constant, there will be no action taken on this species
                        if (// not a constant
                        !simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
                            int stoi = ((Product) reacPart[j]).getStoichiometry();
                            action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-stoi));
                            jp.addAction(action);
                        }
                    }
                }
                // add jump process to compartment subDomain
                subDomain.addJumpProcess(jp);
            }
        // end of if(isForwardRateNonZero), if(isReverseRateNonRate)
        } else if (// flux reactions
        reactionStep instanceof FluxReaction) {
            // we could set jump processes for general flux rate in forms of p1*Sout + p2*Sin
            if (kinetics.getKineticsDescription().equals(KineticsDescription.General) || kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
                Expression fluxRate = new Expression(kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate), reactionStep.getNameScope());
                // we have to pass the math description para to flux solver, coz somehow math description in simulation context is not updated.
                // forward and reverse rate parameters may be null
                Parameter forwardRateParameter = null;
                Parameter reverseRateParameter = null;
                if (kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
                    forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
                    reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
                }
                MassActionSolver.MassActionFunction fluxFunc = MassActionSolver.solveMassAction(forwardRateParameter, reverseRateParameter, fluxRate, (FluxReaction) reactionStep);
                // create jump process for forward flux if it exists.
                Expression rsStructureSize = new Expression(reactionStep.getStructure().getStructureSize(), getNameScope());
                VCUnitDefinition probRateUnit = modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getAreaUnit()).divideBy(modelUnitSystem.getTimeUnit());
                Expression rsRateUnitFactor = getUnitFactor(probRateUnit.divideBy(modelUnitSystem.getFluxReactionUnit()));
                if (fluxFunc.getForwardRate() != null && !fluxFunc.getForwardRate().isZero()) {
                    Expression rate = fluxFunc.getForwardRate();
                    // get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
                    if (fluxFunc.getReactants().size() != 1) {
                        throw new MappingException("Flux " + reactionStep.getName() + " should have only one reactant.");
                    }
                    SpeciesContext scReactant = fluxFunc.getReactants().get(0).getSpeciesContext();
                    Expression scConcExpr = new Expression(getSpeciesConcentrationParameter(scReactant), getNameScope());
                    Expression probExp = Expression.mult(rate, rsRateUnitFactor, rsStructureSize, scConcExpr);
                    // jump process name
                    // +"_reverse";
                    String jpName = TokenMangler.mangleToSName(reactionStep.getName());
                    ProbabilityParameter probParm = null;
                    try {
                        probParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, probExp, PARAMETER_ROLE_P, probabilityParamUnit, reactionStep);
                    } catch (PropertyVetoException pve) {
                        pve.printStackTrace();
                        throw new MappingException(pve.getMessage());
                    }
                    // add probability to function or constant
                    String ms = getMathSymbol(probParm, geometryClass);
                    Expression is = getIdentifierSubstitutions(probExp, probabilityParamUnit, geometryClass);
                    Variable nfoc = newFunctionOrConstant(ms, is, geometryClass);
                    varHash.addVariable(nfoc);
                    JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, geometryClass)));
                    // actions
                    Action action = null;
                    SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-1));
                        jp.addAction(action);
                    }
                    sc = fluxFunc.getProducts().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(1));
                        jp.addAction(action);
                    }
                    subDomain.addJumpProcess(jp);
                }
                // create jump process for reverse flux if it exists.
                if (fluxFunc.getReverseRate() != null && !fluxFunc.getReverseRate().isZero()) {
                    // jump process name
                    String jpName = TokenMangler.mangleToSName(reactionStep.getName()) + PARAMETER_PROBABILITY_RATE_REVERSE_SUFFIX;
                    Expression rate = fluxFunc.getReverseRate();
                    // get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
                    if (fluxFunc.getProducts().size() != 1) {
                        throw new MappingException("Flux " + reactionStep.getName() + " should have only one product.");
                    }
                    SpeciesContext scProduct = fluxFunc.getProducts().get(0).getSpeciesContext();
                    Expression scConcExpr = new Expression(getSpeciesConcentrationParameter(scProduct), getNameScope());
                    Expression probExp = Expression.mult(rate, rsRateUnitFactor, rsStructureSize, scConcExpr);
                    ProbabilityParameter probRevParm = null;
                    try {
                        probRevParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, probExp, PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionStep);
                    } catch (PropertyVetoException pve) {
                        pve.printStackTrace();
                        throw new MappingException(pve.getMessage());
                    }
                    // add probability to function or constant
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, geometryClass), getIdentifierSubstitutions(probExp, probabilityParamUnit, geometryClass), geometryClass));
                    JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, geometryClass)));
                    // actions
                    Action action = null;
                    SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(1));
                        jp.addAction(action);
                    }
                    sc = fluxFunc.getProducts().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-1));
                        jp.addAction(action);
                    }
                    subDomain.addJumpProcess(jp);
                }
            }
        }
    // end of if (simplereaction)...else if(fluxreaction)
    }
// end of reaction step loop
}

Example 5 with SimpleReaction

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

the class StochMathMapping_4_8 method refreshSpeciesContextMappings.

/**
 * Insert the method's description here.
 * Creation date: (10/26/2006 11:47:26 AM)
 * @exception cbit.vcell.parser.ExpressionException The exception description.
 * @exception cbit.vcell.mapping.MappingException The exception description.
 * @exception cbit.vcell.math.MathException The exception description.
 */
private void refreshSpeciesContextMappings() throws cbit.vcell.parser.ExpressionException, MappingException, cbit.vcell.math.MathException {
    // 
    // create a SpeciesContextMapping for each speciesContextSpec.
    // 
    // set initialExpression from SpeciesContextSpec.
    // set diffusing5
    // set variable (only if "Constant" or "Function", else leave it as null)-----why commented?
    // 
    // 
    // have to put geometric paras into mathsymbolmapping, since species initial condition needs the volume size symbol.
    // and the parameters later on were added into contants or functions in refreshMathDescription()
    // 
    StructureMapping[] structureMappings = getSimulationContext().getGeometryContext().getStructureMappings();
    for (int i = 0; i < structureMappings.length; i++) {
        StructureMapping sm = structureMappings[i];
        StructureMapping.StructureMappingParameter parm = sm.getParameterFromRole(StructureMapping.ROLE_Size);
        getMathSymbol(parm, sm);
    }
    getSpeciesContextMappingList().removeAllElements();
    SpeciesContextSpec[] speciesContextSpecs = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        SpeciesContextSpec scs = speciesContextSpecs[i];
        SpeciesContextMapping scm = new SpeciesContextMapping(scs.getSpeciesContext());
        scm.setPDERequired(getSimulationContext().isPDERequired(scs.getSpeciesContext()));
        scm.setHasEventAssignment(getSimulationContext().hasEventAssignment(scs.getSpeciesContext()));
        scm.setHasHybridReaction(false);
        for (ReactionSpec reactionSpec : getSimulationContext().getReactionContext().getReactionSpecs()) {
            if (!reactionSpec.isExcluded() && reactionSpec.hasHybrid(getSimulationContext(), scs.getSpeciesContext())) {
                scm.setHasHybridReaction(true);
            }
        }
        // scm.setAdvecting(isAdvectionRequired(scs.getSpeciesContext()));
        if (scs.isConstant()) {
            SpeciesContextSpec.SpeciesContextSpecParameter initCountParm = scs.getInitialCountParameter();
            SpeciesContextSpec.SpeciesContextSpecParameter initConcParm = scs.getInitialConcentrationParameter();
            Expression initCondInCount = null;
            // initial condition is concentration
            if (initConcParm != null && initConcParm.getExpression() != null) {
                initCondInCount = getExpressionConcToAmt(new Expression(initConcParm, getNameScope()), speciesContextSpecs[i].getSpeciesContext());
            } else {
                initCondInCount = new Expression(initCountParm, getNameScope());
            }
            // initCondInCount.bindExpression(this);
            initCondInCount = getSubstitutedExpr(initCondInCount, true, true);
            scm.setDependencyExpression(initCondInCount);
        }
        // 
        // test if participant in fast reaction step, request elimination if possible
        // 
        scm.setFastParticipant(false);
        ReactionSpec[] reactionSpecs = getSimulationContext().getReactionContext().getReactionSpecs();
        for (int j = 0; j < reactionSpecs.length; j++) {
            ReactionSpec reactionSpec = reactionSpecs[j];
            if (reactionSpec.isExcluded()) {
                continue;
            }
            ReactionStep rs = reactionSpec.getReactionStep();
            if (rs instanceof SimpleReaction && rs.countNumReactionParticipants(scs.getSpeciesContext()) > 0) {
                if (reactionSpec.isFast()) {
                    scm.setFastParticipant(true);
                }
            }
        }
        getSpeciesContextMappingList().addElement(scm);
    }
}