Search in sources :

Example 6 with MacroscopicRateConstant

use of cbit.vcell.math.MacroscopicRateConstant in project vcell by virtualcell.

the class RulebasedMathMapping method addStrictMassActionParticleJumpProcess.

private void addStrictMassActionParticleJumpProcess(VariableHash varHash, GeometryClass geometryClass, SubDomain subDomain, ReactionRule reactionRule, String jpName, ArrayList<ParticleVariable> reactantParticles, ArrayList<ParticleVariable> productParticles, ArrayList<Action> forwardActions, ArrayList<Action> reverseActions) throws ExpressionException, ExpressionBindingException, PropertyVetoException, MathException, MappingException {
    String reactionRuleName = reactionRule.getName();
    RbmKineticLaw kinetics = reactionRule.getKineticLaw();
    RulebasedTransformation ruleBasedTransformation = ((RulebasedTransformation) getTransformation());
    if (kinetics.getRateLawType() != RbmKineticLaw.RateLawType.MassAction) {
        throw new RuntimeException("expecting mass action kinetics for reaction rule " + reactionRuleName);
    }
    // 
    // construct stochastic forward or reverse rate expression (separately).  Transform from
    // original expression of "concentrationRate" in terms of rateParameter and reactants/products in concentrations
    // to
    // new stochastic expression of "molecularRate" in terms of forwardRateParameter, reactants/products in molecules, structure sizes, and unit conversions.
    // 
    // (1)  concentrationRate = K * [s0] * [s1]    [uM.s-1]  or   [molecules.um-3.s-1]   or   [molecules.um-2.s-1]  (or other)
    // (2)  molecularRate = P * <s0> * <s1>        [molecules.s-1]
    // 
    // in this math description, we are using <s_i> [molecules], but original kinetics were in [s_i] [uM or molecules.um-2].
    // so through a change in variable to get things in terms of <s_i>.  <<<< Here P is the desired stochastic rate coefficient. >>>
    // 
    // (3)  let [s_i] = <s_i>/structsize(s_i)*unitConversionFactor(substanceunit([s_i])/substanceunit(<s_i>))
    // 
    // in addition to the change in variables, we need to transform the entire expression from concentration/time to molecules/time
    // 
    // (4)  let molecularRate = concentrationRate * structSize(reaction) * unitConversionFactor(substanceunit(molecularRate)/substanceunit(concentrationRate))
    // 
    // (5)  in general, concentationRate = K * PRODUCT([s_i])
    // 
    // change of variables into stochastic variables used in MathDescription, substituting (3) into (5)
    // 
    // (6)  concentrationRate = K * PRODUCT(<s_i>/structsize(s_i)*unitConversionFactor(substanceunit([s_i])/substanceunit(<s_i>)))
    // 
    // reordering to separate the sizes, the unit conversions and the <s_i>
    // 
    // (7)  concentrationRate = K * PRODUCT(<s_i>) * PRODUCT(1/structsize(s_i)) * unitConversionFactor(PRODUCT(substanceunit([s_i])/substanceunit(<s_i>)))
    // 
    // combining (4) and (7)
    // 
    // (8) molecularRate = K * PRODUCT(<s_i>) * PRODUCT(1/structsize(s_i)) * unitConversionFactor(PRODUCT(substanceunit([s_i])/substanceunit(<s_i>))) * structSize(reaction) * unitConversionFactor(substanceunit(molecularRate)/substanceunit(concentrationRate))
    // 
    // collecting terms of sizes and unit conversions
    // 
    // (9)  molecularRate = K * PRODUCT(<s_i>) * structSize(reaction) / PRODUCT(structsize(s_i)) * unitConversionFactor(substanceunit(molecularRate)/substanceunit(concentrationRate) * PRODUCT(substanceunit([s_i])/substanceunit(<s_i>)))
    // 
    // (10) molecularRate = K * PRODUCT(<s_i>) * sizeFactor * unitConversionFactor(substanceConversionUnit)
    // 
    // where
    // 
    // (11) sizeFactor = structSize(reaction) / PRODUCT(structsize(s_i))
    // (12) substanceConversionUnit = substanceunit(molecularRate)/substanceunit(concentrationRate) * PRODUCT(substanceunit([s_i])/substanceunit(<s_i>))
    // 
    // The ParticleJumpCondition wants a single new rate stochastic, P from equation (2).  Note that PRODUCT(<s_i>) will be captured separately the the reactantPatterns.
    // comparing (2) and (10) we have found P.
    // 
    // (13) P = K * sizeFactor * unitConversionFactor(substanceConversionUnit)
    // 
    // the framework also needs the proper unit for P
    // 
    // (14) Unit(P) = Unit(K) * Unit(sizeFactor) * substanceConversionUnit
    // 
    // 
    ModelUnitSystem modelUnitSystem = getSimulationContext().getModel().getUnitSystem();
    VCUnitDefinition stochasticSubstanceUnit = modelUnitSystem.getStochasticSubstanceUnit();
    VCUnitDefinition reactionRuleSubstanceUnit = modelUnitSystem.getSubstanceUnit(reactionRule.getStructure());
    int forwardRuleIndex = 0;
    // 
    // get forward rate parameter and make sure it is constant valued.
    // 
    Parameter forward_rateParameter = kinetics.getLocalParameter(RbmKineticLawParameterType.MassActionForwardRate);
    Expression substitutedForwardRate = MathUtilities.substituteModelParameters(forward_rateParameter.getExpression(), reactionRule.getNameScope().getScopedSymbolTable());
    if (!substitutedForwardRate.flatten().isNumeric()) {
        throw new MappingException("forward rate constant for reaction rule " + reactionRule.getName() + " is not constant");
    }
    // 
    // create forward sizeExp and forward unitFactor
    // 
    VCUnitDefinition forward_substanceConversionUnit = stochasticSubstanceUnit.divideBy(reactionRuleSubstanceUnit);
    VCUnitDefinition forward_sizeFactorUnit = reactionRule.getStructure().getStructureSize().getUnitDefinition();
    Expression forward_sizeFactor = new Expression(reactionRule.getStructure().getStructureSize(), getNameScope());
    for (ReactantPattern reactantPattern : reactionRule.getReactantPatterns()) {
        Expression reactantSizeExp = new Expression(reactantPattern.getStructure().getStructureSize(), getNameScope());
        VCUnitDefinition reactantSizeUnit = reactantPattern.getStructure().getStructureSize().getUnitDefinition();
        VCUnitDefinition reactantSubstanceUnit = modelUnitSystem.getSubstanceUnit(reactantPattern.getStructure());
        forward_sizeFactor = Expression.div(forward_sizeFactor, reactantSizeExp);
        forward_sizeFactorUnit = forward_sizeFactorUnit.divideBy(reactantSizeUnit);
        forward_substanceConversionUnit = forward_substanceConversionUnit.multiplyBy(reactantSubstanceUnit).divideBy(stochasticSubstanceUnit);
    }
    // simplify sizeFactor (often has size/size/size)
    try {
        forward_sizeFactor = RationalExpUtils.getRationalExp(forward_sizeFactor).simplifyAsExpression();
        forward_sizeFactor.bindExpression(getSimulationContext().getModel());
    } catch (ParseException e) {
        e.printStackTrace();
    }
    Expression forward_rateExp = Expression.mult(new Expression(forward_rateParameter, getNameScope()), forward_sizeFactor, getUnitFactor(forward_substanceConversionUnit)).flatten();
    VCUnitDefinition forward_rateUnit = forward_rateParameter.getUnitDefinition().multiplyBy(forward_sizeFactorUnit).multiplyBy(forward_substanceConversionUnit);
    ProbabilityParameter forward_probParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, forward_rateExp, PARAMETER_ROLE_P, forward_rateUnit, reactionRule);
    // add probability to function or constant
    varHash.addVariable(newFunctionOrConstant(getMathSymbol(forward_probParm, geometryClass), getIdentifierSubstitutions(forward_rateExp, forward_rateUnit, geometryClass), geometryClass));
    // add forward ParticleJumpProcess
    String forward_name = reactionRuleName;
    Expression forward_rate = getIdentifierSubstitutions(new Expression(forward_probParm, getNameScope()), forward_probParm.getUnitDefinition(), geometryClass);
    JumpProcessRateDefinition forward_rateDefinition = new MacroscopicRateConstant(forward_rate);
    ReactionRuleAnalysisReport rrarBiomodelForward = ruleBasedTransformation.getRulesForwardMap().get(reactionRule);
    ProcessSymmetryFactor forwardSymmetryFactor = new ProcessSymmetryFactor(rrarBiomodelForward.getSymmetryFactor());
    ParticleJumpProcess forward_particleJumpProcess = new ParticleJumpProcess(forward_name, reactantParticles, forward_rateDefinition, forwardActions, forwardSymmetryFactor);
    subDomain.addParticleJumpProcess(forward_particleJumpProcess);
    // 
    for (ReactionRule rr : getSimulationContext().getModel().getRbmModelContainer().getReactionRuleList()) {
        if (rr == reactionRule) {
            break;
        }
        forwardRuleIndex++;
        if (rr.isReversible()) {
            forwardRuleIndex++;
        }
    }
    // 
    if (reactionRule.isReversible()) {
        Parameter reverse_rateParameter = kinetics.getLocalParameter(RbmKineticLawParameterType.MassActionReverseRate);
        if (reverse_rateParameter == null || reverse_rateParameter.getExpression() == null) {
            throw new MappingException("reverse rate constant for reaction rule " + reactionRule.getName() + " is missing");
        }
        {
            Expression substitutedReverseRate = MathUtilities.substituteModelParameters(reverse_rateParameter.getExpression(), reactionRule.getNameScope().getScopedSymbolTable());
            if (!substitutedReverseRate.flatten().isNumeric()) {
                throw new MappingException("reverse rate constant for reaction rule " + reactionRule.getName() + " is not constant");
            }
        }
        // 
        // create reverse sizeExp and reverse unitFactor
        // 
        VCUnitDefinition reverse_substanceConversionUnit = stochasticSubstanceUnit.divideBy(reactionRuleSubstanceUnit);
        VCUnitDefinition reverse_sizeFactorUnit = reactionRule.getStructure().getStructureSize().getUnitDefinition();
        Expression reverse_sizeFactor = new Expression(reactionRule.getStructure().getStructureSize(), getNameScope());
        for (ProductPattern productPattern : reactionRule.getProductPatterns()) {
            Expression reactantSizeExp = new Expression(productPattern.getStructure().getStructureSize(), getNameScope());
            VCUnitDefinition reactantSizeUnit = productPattern.getStructure().getStructureSize().getUnitDefinition();
            VCUnitDefinition reactantSubstanceUnit = modelUnitSystem.getSubstanceUnit(productPattern.getStructure());
            reverse_sizeFactor = Expression.div(reverse_sizeFactor, reactantSizeExp);
            reverse_sizeFactorUnit = reverse_sizeFactorUnit.divideBy(reactantSizeUnit);
            reverse_substanceConversionUnit = reverse_substanceConversionUnit.multiplyBy(reactantSubstanceUnit).divideBy(stochasticSubstanceUnit);
        }
        // simplify sizeFactor (often has size/size/size)
        try {
            reverse_sizeFactor = RationalExpUtils.getRationalExp(reverse_sizeFactor).simplifyAsExpression();
            reverse_sizeFactor.bindExpression(getSimulationContext().getModel());
        } catch (ParseException e) {
            e.printStackTrace();
        }
        Expression reverse_rateExp = Expression.mult(new Expression(reverse_rateParameter, getNameScope()), reverse_sizeFactor, getUnitFactor(reverse_substanceConversionUnit)).flatten();
        VCUnitDefinition reverse_rateUnit = reverse_rateParameter.getUnitDefinition().multiplyBy(reverse_sizeFactorUnit).multiplyBy(reverse_substanceConversionUnit);
        // if the reaction has forward rate (Mass action,HMMs), or don't have either forward or reverse rate (some other rate laws--like general)
        // we process it as forward reaction
        // get jump process name
        ProbabilityParameter reverse_probParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName + "_reverse", reverse_rateExp, PARAMETER_ROLE_P_reverse, reverse_rateUnit, reactionRule);
        // add probability to function or constant
        varHash.addVariable(newFunctionOrConstant(getMathSymbol(reverse_probParm, geometryClass), getIdentifierSubstitutions(reverse_rateExp, reverse_rateUnit, geometryClass), geometryClass));
        // add reverse ParticleJumpProcess
        Expression reverse_rate = getIdentifierSubstitutions(new Expression(reverse_probParm, getNameScope()), reverse_probParm.getUnitDefinition(), geometryClass);
        String reverse_name = reactionRuleName + "_reverse";
        JumpProcessRateDefinition reverse_rateDefinition = new MacroscopicRateConstant(reverse_rate);
        ReactionRuleAnalysisReport rrarBiomodelReverse = ruleBasedTransformation.getRulesReverseMap().get(reactionRule);
        ProcessSymmetryFactor reverseSymmetryFactor = new ProcessSymmetryFactor(rrarBiomodelReverse.getSymmetryFactor());
        ParticleJumpProcess reverse_particleJumpProcess = new ParticleJumpProcess(reverse_name, productParticles, reverse_rateDefinition, reverseActions, reverseSymmetryFactor);
        subDomain.addParticleJumpProcess(reverse_particleJumpProcess);
        // 
        // check reverse direction mapping and operations with RuleAnalysis.
        // 
        int reverseRuleIndex = forwardRuleIndex + 1;
        ReactionRuleAnalysisReport rrar = ruleBasedTransformation.getRulesReverseMap().get(reactionRule);
        jumpProcessMap.put(reverse_particleJumpProcess, rrar);
    }
}
Also used : ReactionRuleAnalysisReport(cbit.vcell.mapping.RulebasedTransformer.ReactionRuleAnalysisReport) JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition) ReactionRule(cbit.vcell.model.ReactionRule) ProductPattern(cbit.vcell.model.ProductPattern) ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess) RbmKineticLaw(cbit.vcell.model.RbmKineticLaw) ProcessSymmetryFactor(cbit.vcell.math.ParticleJumpProcess.ProcessSymmetryFactor) VCUnitDefinition(cbit.vcell.units.VCUnitDefinition) Expression(cbit.vcell.parser.Expression) Parameter(cbit.vcell.model.Parameter) UnresolvedParameter(cbit.vcell.mapping.ParameterContext.UnresolvedParameter) LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter) ModelParameter(cbit.vcell.model.Model.ModelParameter) MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant) ParseException(jscl.text.ParseException) RulebasedTransformation(cbit.vcell.mapping.RulebasedTransformer.RulebasedTransformation) ModelUnitSystem(cbit.vcell.model.ModelUnitSystem) ReactantPattern(cbit.vcell.model.ReactantPattern)

Example 7 with MacroscopicRateConstant

use of cbit.vcell.math.MacroscopicRateConstant in project vcell by virtualcell.

the class Xmlproducer method getXML.

private org.jdom.Element getXML(ParticleJumpProcess param) {
    org.jdom.Element particleJumpProcessElement = new org.jdom.Element(XMLTags.ParticleJumpProcessTag);
    // name
    particleJumpProcessElement.setAttribute(XMLTags.NameAttrTag, mangle(param.getName()));
    if (param.getProcessSymmetryFactor() != null) {
        particleJumpProcessElement.setAttribute(XMLTags.ProcessSymmetryFactorAttrTag, Double.toString(param.getProcessSymmetryFactor().getFactor()));
    }
    // Selected Particle
    for (ParticleVariable vpv : param.getParticleVariables()) {
        Element e = new Element(XMLTags.SelectedParticleTag);
        e.setAttribute(XMLTags.NameAttrTag, mangle(vpv.getName()));
        particleJumpProcessElement.addContent(e);
    }
    // probability rate
    Element prob = null;
    JumpProcessRateDefinition particleProbabilityRate = param.getParticleRateDefinition();
    if (particleProbabilityRate instanceof MacroscopicRateConstant) {
        prob = new Element(XMLTags.MacroscopicRateConstantTag);
        prob.addContent(mangleExpression(((MacroscopicRateConstant) particleProbabilityRate).getExpression()));
    } else if (particleProbabilityRate instanceof InteractionRadius) {
        prob = new Element(XMLTags.InteractionRadiusTag);
        prob.addContent(mangleExpression(((InteractionRadius) particleProbabilityRate).getExpression()));
    } else {
        throw new RuntimeException("ParticleRateDefinition in XmlProducer not implemented");
    }
    particleJumpProcessElement.addContent(prob);
    // Actions
    for (Action action : param.getActions()) {
        particleJumpProcessElement.addContent(getXML(action));
    }
    return particleJumpProcessElement;
}
Also used : JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition) Action(cbit.vcell.math.Action) InteractionRadius(cbit.vcell.math.InteractionRadius) VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable) MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable) ParticleVariable(cbit.vcell.math.ParticleVariable) Element(org.jdom.Element) MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant) Element(org.jdom.Element)

Aggregations

JumpProcessRateDefinition (cbit.vcell.math.JumpProcessRateDefinition)7 MacroscopicRateConstant (cbit.vcell.math.MacroscopicRateConstant)7 Action (cbit.vcell.math.Action)6 ParticleJumpProcess (cbit.vcell.math.ParticleJumpProcess)6 Expression (cbit.vcell.parser.Expression)6 InteractionRadius (cbit.vcell.math.InteractionRadius)5 MembraneParticleVariable (cbit.vcell.math.MembraneParticleVariable)5 ParticleVariable (cbit.vcell.math.ParticleVariable)5 Variable (cbit.vcell.math.Variable)5 VolumeParticleVariable (cbit.vcell.math.VolumeParticleVariable)5 ArrayList (java.util.ArrayList)5 CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)4 MathException (cbit.vcell.math.MathException)4 ExpressionException (cbit.vcell.parser.ExpressionException)4 MathDescription (cbit.vcell.math.MathDescription)3 MembraneSubDomain (cbit.vcell.math.MembraneSubDomain)3 ProcessSymmetryFactor (cbit.vcell.math.ParticleJumpProcess.ProcessSymmetryFactor)3 SubDomain (cbit.vcell.math.SubDomain)3 Element (org.jdom.Element)3 GeometryClass (cbit.vcell.geometry.GeometryClass)2