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

Example 1 with FluxReaction

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

the class PathwayMapping method createReactionStepsFromTableRow.

private void createReactionStepsFromTableRow(BioModel bioModel, Transport bioPaxObject, double stoich, String id, String location, ArrayList<ConversionTableRow> conversionTableRows, boolean addSubunits) throws Exception {
    if (bioModel == null) {
        return;
    }
    for (Process process : BioPAXUtil.getAllProcesses(bioModel.getPathwayModel(), bioPaxObject)) {
        // use user defined id as the name of the reaction name
        // get participants from table rows
        ArrayList<ConversionTableRow> participants = new ArrayList<ConversionTableRow>();
        for (ConversionTableRow ctr : conversionTableRows) {
            if (ctr.interactionId().equals(bioPaxObject.getID())) {
                participants.add(ctr);
            }
        }
        // create reaction object
        String name = getSafetyName(process.getName() + "_" + location);
        if (bioModel.getModel().getReactionStep(name) == null) {
            // create a new reactionStep object
            FluxReaction fluxReactionStep = bioModel.getModel().createFluxReaction((Membrane) bioModel.getModel().getStructure(location));
            fluxReactionStep.setName(name);
            RelationshipObject newRelationship = new RelationshipObject(fluxReactionStep, bioPaxObject);
            bioModel.getRelationshipModel().addRelationshipObject(newRelationship);
            createReactionStep(bioModel, process, fluxReactionStep, newRelationship, participants, addSubunits);
        } else {
            // bioModel.getModel().getReactionStep(safeId).setStructure(bioModel.getModel().getStructure(location));
            // add missing parts for the existing reactionStep
            RelationshipObject newRelationship = new RelationshipObject(bioModel.getModel().getReactionStep(name), bioPaxObject);
            bioModel.getRelationshipModel().addRelationshipObject(newRelationship);
            createReactionStep(bioModel, process, bioModel.getModel().getReactionStep(name), newRelationship, participants, addSubunits);
        }
    }
}

Example 2 with FluxReaction

use of cbit.vcell.model.FluxReaction 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 3 with FluxReaction

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

the class XmlReader method getFluxReaction.

/**
 * This method returns a FluxReaction object from a XML element.
 * Creation date: (3/16/2001 11:52:02 AM)
 * @return cbit.vcell.model.FluxReaction
 * @param param org.jdom.Element
 * @throws XmlParseException
 * @throws PropertyVetoException
 * @throws ModelException
 * @throws Exception
 */
private FluxReaction getFluxReaction(Element param, Model model) throws XmlParseException, PropertyVetoException {
    // retrieve the key if there is one
    KeyValue key = null;
    String keystring = param.getAttributeValue(XMLTags.KeyValueAttrTag);
    if (keystring != null && keystring.length() > 0 && this.readKeysFlag) {
        key = new KeyValue(keystring);
    }
    // resolve reference to the Membrane
    String structureName = unMangle(param.getAttributeValue(XMLTags.StructureAttrTag));
    Membrane structureref = (Membrane) model.getStructure(structureName);
    if (structureref == null) {
        throw new XmlParseException("The membrane " + structureName + " could not be resolved in the dictionnary!");
    }
    // -- Instantiate new FluxReaction --
    FluxReaction fluxreaction = null;
    String name = unMangle(param.getAttributeValue(XMLTags.NameAttrTag));
    String reversibleAttributeValue = param.getAttributeValue(XMLTags.ReversibleAttrTag);
    boolean bReversible = true;
    if (reversibleAttributeValue != null) {
        if (Boolean.TRUE.toString().equals(reversibleAttributeValue)) {
            bReversible = true;
        } else if (Boolean.FALSE.toString().equals(reversibleAttributeValue)) {
            bReversible = false;
        } else {
            throw new RuntimeException("unexpected value " + reversibleAttributeValue + " for reversible flag for reaction " + name);
        }
    }
    try {
        fluxreaction = new FluxReaction(model, structureref, key, name, bReversible);
        fluxreaction.setModel(model);
    } catch (Exception e) {
        e.printStackTrace();
        throw new XmlParseException("An exception occurred while trying to create the FluxReaction " + name, e);
    }
    // resolve reference to the fluxCarrier
    if (param.getAttribute(XMLTags.FluxCarrierAttrTag) != null) {
        String speciesname = unMangle(param.getAttributeValue(XMLTags.FluxCarrierAttrTag));
        Species specieref = model.getSpecies(speciesname);
        if (specieref != null) {
            Feature insideFeature = model.getStructureTopology().getInsideFeature(structureref);
            try {
                if (insideFeature != null) {
                    SpeciesContext insideSpeciesContext = model.getSpeciesContext(specieref, insideFeature);
                    fluxreaction.addProduct(insideSpeciesContext, 1);
                }
                Feature outsideFeature = model.getStructureTopology().getOutsideFeature(structureref);
                if (outsideFeature != null) {
                    SpeciesContext outsideSpeciesContext = model.getSpeciesContext(specieref, outsideFeature);
                    fluxreaction.addReactant(outsideSpeciesContext, 1);
                }
            } catch (ModelException e) {
                e.printStackTrace(System.out);
                throw new XmlParseException(e.getMessage());
            }
        }
    }
    // Annotation
    // String rsAnnotation = null;
    // String annotationText = param.getChildText(XMLTags.AnnotationTag, vcNamespace);
    // if (annotationText!=null && annotationText.length()>0) {
    // rsAnnotation = unMangle(annotationText);
    // }
    // fluxreaction.setAnnotation(rsAnnotation);
    // set the fluxOption
    String fluxOptionString = null;
    fluxOptionString = param.getAttributeValue(XMLTags.FluxOptionAttrTag);
    if (fluxOptionString != null && fluxOptionString.length() > 0) {
        try {
            if (fluxOptionString.equals(XMLTags.FluxOptionElectricalOnly)) {
                fluxreaction.setPhysicsOptions(FluxReaction.PHYSICS_ELECTRICAL_ONLY);
            } else if (fluxOptionString.equals(XMLTags.FluxOptionMolecularAndElectrical)) {
                fluxreaction.setPhysicsOptions(FluxReaction.PHYSICS_MOLECULAR_AND_ELECTRICAL);
            } else if (fluxOptionString.equals(XMLTags.FluxOptionMolecularOnly)) {
                fluxreaction.setPhysicsOptions(FluxReaction.PHYSICS_MOLECULAR_ONLY);
            }
        } catch (java.beans.PropertyVetoException e) {
            e.printStackTrace(System.out);
            throw new XmlParseException("A propertyVetoException was fired when setting the fluxOption to the flux reaction " + name, e);
        }
    }
    // Add Reactants, if any
    try {
        Iterator<Element> iterator = param.getChildren(XMLTags.ReactantTag, vcNamespace).iterator();
        while (iterator.hasNext()) {
            Element temp = iterator.next();
            // Add Reactant to this SimpleReaction
            fluxreaction.addReactionParticipant(getReactant(temp, fluxreaction, model));
        }
    } catch (java.beans.PropertyVetoException e) {
        e.printStackTrace();
        throw new XmlParseException("Error adding a reactant to the reaction " + name + " : " + e.getMessage());
    }
    // Add Products, if any
    try {
        Iterator<Element> iterator = param.getChildren(XMLTags.ProductTag, vcNamespace).iterator();
        while (iterator.hasNext()) {
            Element temp = iterator.next();
            // Add Product to this simplereaction
            fluxreaction.addReactionParticipant(getProduct(temp, fluxreaction, model));
        }
    } catch (java.beans.PropertyVetoException e) {
        e.printStackTrace();
        throw new XmlParseException("Error adding a product to the reaction " + name + " : " + e.getMessage());
    }
    // Add Catalyst(Modifiers) (if there are)
    Iterator<Element> iterator = param.getChildren(XMLTags.CatalystTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element temp = iterator.next();
        fluxreaction.addReactionParticipant(getCatalyst(temp, fluxreaction, model));
    }
    // Add Kinetics
    fluxreaction.setKinetics(getKinetics(param.getChild(XMLTags.KineticsTag, vcNamespace), fluxreaction, model));
    // set the valence (for legacy support for "chargeCarrierValence" stored with reaction).
    String valenceString = null;
    try {
        valenceString = unMangle(param.getAttributeValue(XMLTags.FluxCarrierValenceAttrTag));
        if (valenceString != null && valenceString.length() > 0) {
            KineticsParameter chargeValenceParameter = fluxreaction.getKinetics().getChargeValenceParameter();
            if (chargeValenceParameter != null) {
                chargeValenceParameter.setExpression(new Expression(Integer.parseInt(unMangle(valenceString))));
            }
        }
    } catch (NumberFormatException e) {
        e.printStackTrace();
        throw new XmlParseException("A NumberFormatException was fired when setting the (integer) valence '" + valenceString + "' (integer) to the flux reaction " + name, e);
    }
    return fluxreaction;
}

Example 4 with FluxReaction

use of cbit.vcell.model.FluxReaction 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 FluxReaction

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

the class StructureAnalyzer method getReactionRateExpression.

public Expression getReactionRateExpression(ReactionStep reactionStep, ReactionParticipant reactionParticipant) throws Exception {
    if (reactionParticipant instanceof Catalyst) {
        throw new Exception("Catalyst " + reactionParticipant + " doesn't have a rate for this reaction");
    // return new Expression(0.0);
    }
    double stoich = reactionStep.getStoichiometry(reactionParticipant.getSpeciesContext());
    if (stoich == 0.0) {
        return new Expression(0.0);
    }
    if (reactionStep.getKinetics() instanceof DistributedKinetics) {
        DistributedKinetics distributedKinetics = (DistributedKinetics) reactionStep.getKinetics();
        if (stoich != 1) {
            Expression exp = Expression.mult(new Expression(stoich), new Expression(distributedKinetics.getReactionRateParameter(), mathMapping_4_8.getNameScope()));
            return exp;
        } else {
            Expression exp = new Expression(distributedKinetics.getReactionRateParameter(), mathMapping_4_8.getNameScope());
            return exp;
        }
    } else if (reactionStep.getKinetics() instanceof LumpedKinetics) {
        Structure.StructureSize structureSize = reactionStep.getStructure().getStructureSize();
        // 
        // need to put this into concentration/time with respect to structure for reaction.
        // 
        LumpedKinetics lumpedKinetics = (LumpedKinetics) reactionStep.getKinetics();
        Expression factor = null;
        ModelUnitSystem unitSystem = mathMapping_4_8.getSimulationContext().getModel().getUnitSystem();
        if (reactionStep.getStructure() instanceof Feature || ((reactionStep.getStructure() instanceof Membrane) && reactionStep instanceof FluxReaction)) {
            VCUnitDefinition lumpedToVolumeSubstance = unitSystem.getVolumeSubstanceUnit().divideBy(unitSystem.getLumpedReactionSubstanceUnit());
            factor = Expression.div(new Expression(lumpedToVolumeSubstance.getDimensionlessScale()), new Expression(structureSize, mathMapping_4_8.getNameScope()));
        } else if (reactionStep.getStructure() instanceof Membrane && reactionStep instanceof SimpleReaction) {
            VCUnitDefinition lumpedToVolumeSubstance = unitSystem.getMembraneSubstanceUnit().divideBy(unitSystem.getLumpedReactionSubstanceUnit());
            factor = Expression.div(new Expression(lumpedToVolumeSubstance.getDimensionlessScale()), new Expression(structureSize, mathMapping_4_8.getNameScope()));
        } else {
            throw new RuntimeException("failed to create reaction rate expression for reaction " + reactionStep.getName() + ", with kinetic type of " + reactionStep.getKinetics().getClass().getName());
        }
        if (stoich != 1) {
            Expression exp = Expression.mult(new Expression(stoich), Expression.mult(new Expression(lumpedKinetics.getLumpedReactionRateParameter(), mathMapping_4_8.getNameScope()), factor));
            return exp;
        } else {
            Expression exp = Expression.mult(new Expression(lumpedKinetics.getLumpedReactionRateParameter(), mathMapping_4_8.getNameScope()), factor);
            return exp;
        }
    } else {
        throw new RuntimeException("unexpected kinetic type " + reactionStep.getKinetics().getClass().getName());
    }
}