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

Example 56 with Membrane

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

the class MathMapping_4_8 method refreshStructureAnalyzers.

/**
 * This method was created by a SmartGuide.
 */
protected void refreshStructureAnalyzers() {
    structureAnalyzerList.removeAllElements();
    // 
    // update structureAnalyzer list if any subVolumes were added
    // 
    SubVolume[] subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
    for (int j = 0; j < subVolumes.length; j++) {
        SubVolume subVolume = (SubVolume) subVolumes[j];
        if (getVolumeStructureAnalyzer(subVolume) == null) {
            structureAnalyzerList.addElement(new VolumeStructureAnalyzer(this, subVolume));
        }
        // 
        // Add a MembraneStructureAnalyzer if necessary
        // 
        // go through list of MembraneMappings and determine if inner and outer compartment
        // are both mapped to subVolumes, then add
        // 
        Structure[] structures = getStructures(subVolume);
        if (structures != null) {
            for (int i = 0; i < structures.length; i++) {
                if (structures[i] instanceof Membrane) {
                    Membrane membrane = (Membrane) structures[i];
                    MembraneMapping mm = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membrane);
                    if (mm != null) {
                        if (getResolved(mm) && getMembraneStructureAnalyzer(membrane) == null) {
                            SubVolume outerSubVolume = getSubVolume(((FeatureMapping) simContext.getGeometryContext().getStructureMapping(simContext.getModel().getStructureTopology().getOutsideFeature(membrane))));
                            structureAnalyzerList.addElement(new MembraneStructureAnalyzer(this, membrane, subVolume, outerSubVolume));
                        }
                    }
                }
            }
        }
    }
    // 
    // invoke all structuralAnalyzers
    // 
    Enumeration<StructureAnalyzer> enum1 = getStructureAnalyzers();
    while (enum1.hasMoreElements()) {
        StructureAnalyzer sa = enum1.nextElement();
        sa.refresh();
    }
}

Example 57 with Membrane

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

the class MathMapping_4_8 method refreshVariables.

/**
 * This method was created in VisualAge.
 */
private void refreshVariables() throws MappingException {
    // System.out.println("MathMapping.refreshVariables()");
    // 
    // non-constant dependent variables require a function
    // 
    Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = enum1.nextElement();
        SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(scm.getSpeciesContext());
        if (scm.getDependencyExpression() != null && !scs.isConstant()) {
            // scm.setVariable(new Function(scm.getSpeciesContext().getName(),scm.getDependencyExpression()));
            scm.setVariable(null);
        }
    }
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = enum1.nextElement();
        SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(scm.getSpeciesContext());
        if (getSimulationContext().hasEventAssignment(scs.getSpeciesContext())) {
            scm.setDependencyExpression(null);
        }
    }
    // 
    // non-constant independent variables require either a membrane or volume variable
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(scm.getSpeciesContext());
        if (scm.getDependencyExpression() == null && (!scs.isConstant() || getSimulationContext().hasEventAssignment(scs.getSpeciesContext()))) {
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
            Structure struct = scm.getSpeciesContext().getStructure();
            if (struct instanceof Feature) {
                if (getResolved(sm)) {
                    scm.setVariable(getResolvedVolVariable(scm.getSpeciesContext().getSpecies()));
                } else {
                    scm.setVariable(new VolVariable(scm.getSpeciesContext().getName(), nullDomain));
                }
            } else if (struct instanceof Membrane) {
                if (getResolved(sm)) {
                    scm.setVariable(new MemVariable(scm.getSpeciesContext().getName(), nullDomain));
                } else {
                    scm.setVariable(new VolVariable(scm.getSpeciesContext().getName(), nullDomain));
                }
            } else {
                throw new MappingException("class " + scm.getSpeciesContext().getStructure().getClass() + " not supported");
            }
            mathSymbolMapping.put(scm.getSpeciesContext(), scm.getVariable().getName());
        }
    }
}

Example 58 with Membrane

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

the class MembraneStructureAnalyzer method refreshResolvedFluxes.

/**
 * This method was created in VisualAge.
 */
void refreshResolvedFluxes() throws Exception {
    // System.out.println("MembraneStructureAnalyzer.refreshResolvedFluxes()");
    GeometryContext geoContext = mathMapping_4_8.getSimulationContext().getGeometryContext();
    StructureTopology structTopology = mathMapping_4_8.getSimulationContext().getModel().getStructureTopology();
    Vector<ResolvedFlux> resolvedFluxList = new Vector<ResolvedFlux>();
    // 
    // for each reaction, get all fluxReactions associated with this membrane
    // 
    Vector<ReactionStep> fluxList = new Vector<ReactionStep>();
    ReactionSpec[] reactionSpecs = mathMapping_4_8.getSimulationContext().getReactionContext().getReactionSpecs();
    for (int j = 0; j < reactionSpecs.length; j++) {
        if (reactionSpecs[j].isExcluded()) {
            continue;
        }
        ReactionStep rs = reactionSpecs[j].getReactionStep();
        if (rs.getStructure() == getMembrane()) {
            if (rs instanceof FluxReaction) {
                fluxList.addElement(rs);
            }
        }
    }
    // 
    for (int i = 0; i < fluxList.size(); i++) {
        FluxReaction fr = (FluxReaction) fluxList.elementAt(i);
        Species fluxCarrier = null;
        for (ReactionParticipant rp : fr.getReactionParticipants()) {
            if (rp instanceof Reactant || rp instanceof Product) {
                if (fluxCarrier == null) {
                    fluxCarrier = rp.getSpecies();
                } else {
                    if (fluxCarrier != rp.getSpecies()) {
                        throw new Exception("Flux reaction '" + fr.getName() + "' with multiple species not allowed in VCell 4.8.");
                    }
                }
            }
        }
        if (fluxCarrier == null) {
            continue;
        }
        ResolvedFlux rf = null;
        for (int j = 0; j < resolvedFluxList.size(); j++) {
            ResolvedFlux rf_tmp = (ResolvedFlux) resolvedFluxList.elementAt(j);
            if (rf_tmp.getSpecies() == fluxCarrier) {
                rf = rf_tmp;
            }
        }
        // 
        // if "inside" speciesContext is not "fixed", add flux to ResolvedFlux
        // 
        SpeciesContext insideSpeciesContext = mathMapping_4_8.getSimulationContext().getModel().getSpeciesContext(fluxCarrier, structTopology.getInsideFeature(getMembrane()));
        SpeciesContextSpec insideSpeciesContextSpec = mathMapping_4_8.getSimulationContext().getReactionContext().getSpeciesContextSpec(insideSpeciesContext);
        // if (!insideSpeciesContextSpec.isConstant()){
        if (bNoFluxIfFixed || !insideSpeciesContextSpec.isConstant()) {
            if (bNoFluxIfFixed && insideSpeciesContextSpec.isConstant()) {
                bNoFluxIfFixedExercised = true;
            }
            if (rf == null) {
                rf = new ResolvedFlux(fluxCarrier, fr.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getUnitDefinition());
                resolvedFluxList.addElement(rf);
            }
            FeatureMapping insideFeatureMapping = (FeatureMapping) geoContext.getStructureMapping((structTopology.getInsideFeature((Membrane) fr.getStructure())));
            Expression residualVolumeFraction = mathMapping_4_8.getResidualVolumeFraction(insideFeatureMapping).renameBoundSymbols(mathMapping_4_8.getNameScope());
            Expression insideFluxCorrection = Expression.invert(residualVolumeFraction);
            // 
            if (bResolvedFluxCorrectionBug && !residualVolumeFraction.compareEqual(new Expression(1.0))) {
                bResolvedFluxCorrectionBugExercised = true;
                System.out.println("MembraneStructureAnalyzer.refreshResolvedFluxes() ... 'ResolvedFluxCorrection' bug compatability mode");
                insideFluxCorrection = new Expression(1.0);
            }
            // 
            if (fr.getKinetics() instanceof DistributedKinetics) {
                Expression reactionRateParameter = new Expression(((DistributedKinetics) fr.getKinetics()).getReactionRateParameter(), mathMapping_4_8.getNameScope());
                if (rf.inFluxExpression.isZero()) {
                    rf.inFluxExpression = Expression.mult(reactionRateParameter, insideFluxCorrection).flatten();
                } else {
                    rf.inFluxExpression = Expression.add(rf.inFluxExpression, Expression.mult(reactionRateParameter, insideFluxCorrection).flatten());
                }
            } else if (fr.getKinetics() instanceof LumpedKinetics) {
                throw new RuntimeException("Lumped Kinetics for fluxes not yet supported");
            } else {
                throw new RuntimeException("unexpected Kinetic type in MembraneStructureAnalyzer.refreshResolvedFluxes()");
            }
        // rf.inFlux.bindExpression(mathMapping);
        }
        SpeciesContext outsideSpeciesContext = mathMapping_4_8.getSimulationContext().getModel().getSpeciesContext(fluxCarrier, structTopology.getOutsideFeature(getMembrane()));
        SpeciesContextSpec outsideSpeciesContextSpec = mathMapping_4_8.getSimulationContext().getReactionContext().getSpeciesContextSpec(outsideSpeciesContext);
        // if (!outsideSpeciesContextSpec.isConstant()){
        if (bNoFluxIfFixed || !outsideSpeciesContextSpec.isConstant()) {
            if (bNoFluxIfFixed && outsideSpeciesContextSpec.isConstant()) {
                bNoFluxIfFixedExercised = true;
            }
            if (rf == null) {
                rf = new ResolvedFlux(fluxCarrier, fr.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getUnitDefinition());
                resolvedFluxList.addElement(rf);
            }
            FeatureMapping outsideFeatureMapping = (FeatureMapping) geoContext.getStructureMapping(structTopology.getOutsideFeature((Membrane) fr.getStructure()));
            Expression residualVolumeFraction = mathMapping_4_8.getResidualVolumeFraction(outsideFeatureMapping).renameBoundSymbols(mathMapping_4_8.getNameScope());
            Expression outsideFluxCorrection = Expression.invert(residualVolumeFraction);
            // 
            if (bResolvedFluxCorrectionBug && !residualVolumeFraction.compareEqual(new Expression(1.0))) {
                bResolvedFluxCorrectionBugExercised = true;
                System.out.println("MembraneStructureAnalyzer.refreshResolvedFluxes() ... 'ResolvedFluxCorrection' bug compatability mode");
                outsideFluxCorrection = new Expression(1.0);
            }
            // 
            if (fr.getKinetics() instanceof DistributedKinetics) {
                Expression reactionRateParameter = new Expression(((DistributedKinetics) fr.getKinetics()).getReactionRateParameter(), mathMapping_4_8.getNameScope());
                if (rf.outFluxExpression.isZero()) {
                    rf.outFluxExpression = Expression.mult(Expression.negate(reactionRateParameter), outsideFluxCorrection).flatten();
                } else {
                    rf.outFluxExpression = Expression.add(rf.outFluxExpression, Expression.mult(Expression.negate(reactionRateParameter), outsideFluxCorrection).flatten());
                }
            } else if (fr.getKinetics() instanceof LumpedKinetics) {
                throw new RuntimeException("Lumped Kinetics not yet supported for Flux Reaction: " + fr.getName());
            } else {
                throw new RuntimeException("unexpected Kinetics type for Flux Reaction " + fr.getName());
            }
        // rf.outFlux.bindExpression(mathMapping);
        }
    }
    // 
    // for each reaction, incorporate all reactionSteps involving binding with volumetric species
    // 
    double kMoleValue = 1 / 602.0;
    for (int i = 0; i < reactionSpecs.length; i++) {
        if (reactionSpecs[i].isExcluded()) {
            continue;
        }
        ReactionStep rs = reactionSpecs[i].getReactionStep();
        if (rs.getStructure() == getMembrane()) {
            if (rs instanceof SimpleReaction) {
                SimpleReaction sr = (SimpleReaction) rs;
                ReactionParticipant[] rp_Array = sr.getReactionParticipants();
                for (int k = 0; k < rp_Array.length; k++) {
                    if (rp_Array[k] instanceof Reactant || rp_Array[k] instanceof Product) {
                        SpeciesContextSpec scs = mathMapping_4_8.getSimulationContext().getReactionContext().getSpeciesContextSpec(rp_Array[k].getSpeciesContext());
                        // if (rp_Array[k].getStructure() instanceof Feature && !scs.isConstant()){
                        if (rp_Array[k].getStructure() instanceof Feature && (bNoFluxIfFixed || !scs.isConstant())) {
                            if (bNoFluxIfFixed && scs.isConstant()) {
                                bNoFluxIfFixedExercised = true;
                            }
                            // 
                            // for each Reactant or Product binding to this membrane...
                            // 
                            // 
                            // get ResolvedFlux for this species
                            // 
                            ResolvedFlux rf = null;
                            for (int j = 0; j < resolvedFluxList.size(); j++) {
                                ResolvedFlux rf_tmp = (ResolvedFlux) resolvedFluxList.elementAt(j);
                                if (rf_tmp.getSpecies() == rp_Array[k].getSpecies()) {
                                    rf = rf_tmp;
                                }
                            }
                            if (rf == null) {
                                rf = new ResolvedFlux(rp_Array[k].getSpecies(), sr.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getUnitDefinition());
                                resolvedFluxList.addElement(rf);
                            }
                            Expression reactionRateExpression = getReactionRateExpression(sr, rp_Array[k]).renameBoundSymbols(mathMapping_4_8.getNameScope());
                            if (rp_Array[k].getStructure() == structTopology.getInsideFeature(getMembrane())) {
                                // 
                                // for binding on inside, add to ResolvedFlux.inFlux
                                // 
                                FeatureMapping insideFeatureMapping = (FeatureMapping) geoContext.getStructureMapping(structTopology.getInsideFeature(getMembrane()));
                                Expression residualVolumeFraction = mathMapping_4_8.getResidualVolumeFraction(insideFeatureMapping).renameBoundSymbols(mathMapping_4_8.getNameScope());
                                Expression insideFluxCorrection = Expression.div(new Expression(kMoleValue), residualVolumeFraction).flatten();
                                // 
                                if (bResolvedFluxCorrectionBug && !residualVolumeFraction.compareEqual(new Expression(1.0))) {
                                    bResolvedFluxCorrectionBugExercised = true;
                                    System.out.println("MembraneStructureAnalyzer.refreshResolvedFluxes() ... 'ResolvedFluxCorrection' bug compatability mode");
                                    insideFluxCorrection = new Expression(kMoleValue);
                                }
                                if (rf.inFluxExpression.isZero()) {
                                    rf.inFluxExpression = Expression.mult(insideFluxCorrection, reactionRateExpression);
                                } else {
                                    rf.inFluxExpression = Expression.add(rf.inFluxExpression, Expression.mult(insideFluxCorrection, reactionRateExpression));
                                }
                            // rf.inFlux.bindExpression(mathMapping);
                            } else if (rp_Array[k].getStructure() == structTopology.getOutsideFeature(getMembrane())) {
                                // 
                                // for binding on outside, add to ResolvedFlux.outFlux
                                // 
                                FeatureMapping outsideFeatureMapping = (FeatureMapping) geoContext.getStructureMapping(structTopology.getOutsideFeature(getMembrane()));
                                Expression residualVolumeFraction = mathMapping_4_8.getResidualVolumeFraction(outsideFeatureMapping).renameBoundSymbols(mathMapping_4_8.getNameScope());
                                Expression outsideFluxCorrection = Expression.div(new Expression(kMoleValue), residualVolumeFraction).flatten();
                                // 
                                if (bResolvedFluxCorrectionBug && !residualVolumeFraction.compareEqual(new Expression(1.0))) {
                                    bResolvedFluxCorrectionBugExercised = true;
                                    System.out.println("MembraneStructureAnalyzer.refreshResolvedFluxes() ... 'ResolvedFluxCorrection' bug compatability mode");
                                    outsideFluxCorrection = new Expression(kMoleValue);
                                }
                                if (rf.outFluxExpression.isZero()) {
                                    rf.outFluxExpression = Expression.mult(outsideFluxCorrection, reactionRateExpression);
                                } else {
                                    rf.outFluxExpression = Expression.add(rf.outFluxExpression, Expression.mult(outsideFluxCorrection, reactionRateExpression));
                                }
                            // rf.outFlux.bindExpression(mathMapping);
                            } else {
                                throw new Exception("SpeciesContext " + rp_Array[k].getSpeciesContext().getName() + " doesn't border membrane " + getMembrane().getName() + " but reacts there");
                            }
                        }
                    }
                }
            }
        }
    }
    // 
    if (resolvedFluxList.size() > 0) {
        resolvedFluxes = new ResolvedFlux[resolvedFluxList.size()];
        resolvedFluxList.copyInto(resolvedFluxes);
    } else {
        resolvedFluxes = null;
    }
}

Example 59 with Membrane

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

the class StochMathMapping_4_8 method getProbabilityRate.

/**
 * Get probability expression for the specific elementary reaction.
 * Input: ReactionStep, the reaction. isForwardDirection, if the elementary reaction is forward from the reactionstep.
 * Output: Expression. the probability expression.
 * Creation date: (9/14/2006 3:22:58 PM)
 */
Expression getProbabilityRate(ReactionStep rs, boolean isForwardDirection) throws MappingException {
    ReactionStep reactionStep = rs;
    Expression probExp = null;
    // get kinetics of the reaction step
    Kinetics kinetics = reactionStep.getKinetics();
    // to compose the rate constant expression e.g. Kf, Kr
    Expression rateConstantExpr = null;
    // to compose the stochastic variable(species) expression, e.g. s*(s-1)*(s-2)* speciesFactor.
    Expression rxnProbabilityExpr = null;
    // to compose the factor that the probability expression multiplies with, which convert the rate expression under stochastic context
    Expression factorExpr = null;
    // the structure where reaction happens
    StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(rs.getStructure());
    Model model = getSimulationContext().getModel();
    try {
        if (// forward reaction
        isForwardDirection) {
            // for HMMs, it's a bit complicated. Vmax/(Km+s)-->Vmax*Size_s/(Km*Size_s+Ns)
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                KineticsParameter kfp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward);
                rateConstantExpr = new Expression(kfp, getNameScope());
            // rateConstantExpr.bindExpression(this);
            }
            // get convert factor for rate constant( membrane:rateConstant*membrane_Size (factor is membrane_size), feature : rateConstant*(feature_size/KMole)(factor is feature_size/KMOLE)) )
            if (sm.getStructure() instanceof Membrane) {
                factorExpr = new Expression(sm.getStructure().getStructureSize(), getNameScope());
            } else {
                factorExpr = new Expression(sm.getStructure().getStructureSize(), getNameScope());
                Expression kmoleExpr = new Expression(1.0 / 602.0);
                factorExpr = Expression.mult(factorExpr, kmoleExpr);
            }
            // complete the probability expression by the reactants' stoichiometries if it is Mass Action rate law
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
                for (int i = 0; i < reacPart.length; i++) {
                    int stoichiometry = 0;
                    if (reacPart[i] instanceof Reactant) {
                        stoichiometry = ((Reactant) reacPart[i]).getStoichiometry();
                        // ******the following part is to form the s*(s-1)(s-2)..(s-stoi+1).portion of the probability rate.
                        StructureMapping reactSM = getSimulationContext().getGeometryContext().getStructureMapping(reacPart[i].getStructure());
                        // factor expression for species
                        Expression speciesFactor = null;
                        // convert speceis' unit from moles/liter to molecules.
                        if (reactSM.getStructure() instanceof Membrane) {
                            speciesFactor = Expression.invert(new Expression(reactSM.getStructure().getStructureSize(), getNameScope()));
                        } else {
                            Expression exp1 = new Expression(1.0 / 602.0);
                            Expression exp2 = new Expression(reactSM.getStructure().getStructureSize(), getNameScope());
                            speciesFactor = Expression.div(Expression.invert(exp1), exp2);
                        }
                        // s*(s-1)(s-2)..(s-stoi+1)
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[i].getSpeciesContext());
                        Expression spCount_exp = new Expression(spCountParam, getNameScope());
                        // species from uM to No. of Particles, form s*(s-1)*(s-2)
                        Expression tempExpr = new Expression(spCount_exp);
                        for (int j = 1; j < stoichiometry; j++) {
                            tempExpr = Expression.mult(tempExpr, Expression.add(spCount_exp, new Expression(-j)));
                        }
                        // update total factor with speceies factor
                        if (stoichiometry == 1) {
                            factorExpr = Expression.mult(factorExpr, speciesFactor);
                        } else if (stoichiometry > 1) {
                            // rxnProbExpr * (structSize^stoichiometry)
                            Expression powerExpr = Expression.power(speciesFactor, new Expression(stoichiometry));
                            factorExpr = Expression.mult(factorExpr, powerExpr);
                        }
                        if (rxnProbabilityExpr == null) {
                            rxnProbabilityExpr = new Expression(tempExpr);
                        } else {
                            // for more than one reactant
                            rxnProbabilityExpr = Expression.mult(rxnProbabilityExpr, tempExpr);
                        }
                    }
                }
            }
        } else // reverse reaction
        {
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                KineticsParameter krp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse);
                rateConstantExpr = new Expression(krp, getNameScope());
            // rateConstantExpr.bindExpression(this);
            }
            // get convert factor for rate constant( membrane:rateConstant*membrane_Size (factor is membrane_size), feature : rateConstant*(feature_size/KMole)(factor is feature_size/KMOLE)) )
            if (sm.getStructure() instanceof Membrane) {
                factorExpr = new Expression(sm.getStructure().getStructureSize(), getNameScope());
            } else {
                factorExpr = new Expression(sm.getStructure().getStructureSize(), getNameScope());
                Expression exp = new Expression(1.0 / 602.0);
                factorExpr = Expression.mult(factorExpr, exp);
            }
            // complete the remaining part of the probability expression by the products' stoichiometries.
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
                for (int i = 0; i < reacPart.length; i++) {
                    int stoichiometry = 0;
                    if (reacPart[i] instanceof Product) {
                        stoichiometry = ((Product) reacPart[i]).getStoichiometry();
                        // ******the following part is to form the s*(s-1)*(s-2)...(s-stoi+1).portion of the probability rate.
                        StructureMapping reactSM = getSimulationContext().getGeometryContext().getStructureMapping(reacPart[i].getStructure());
                        // factor expression for species
                        Expression speciesFactor = null;
                        // convert speceis' unit from moles/liter to molecules.
                        if (reactSM.getStructure() instanceof Membrane) {
                            speciesFactor = Expression.invert(new Expression(reactSM.getStructure().getStructureSize(), getNameScope()));
                        } else {
                            Expression exp1 = new Expression(1.0 / 602.0);
                            Expression exp2 = new Expression(reactSM.getStructure().getStructureSize(), getNameScope());
                            speciesFactor = Expression.div(Expression.invert(exp1), exp2);
                        }
                        // s*(s-1)*(s-2)...(s-stoi+1)
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[i].getSpeciesContext());
                        Expression spCount_exp = new Expression(spCountParam, getNameScope());
                        // species from uM to No. of Particles, form s*(s-1)*(s-2)
                        Expression tempExpr = new Expression(spCount_exp);
                        for (int j = 1; j < stoichiometry; j++) {
                            tempExpr = Expression.mult(tempExpr, Expression.add(spCount_exp, new Expression(-j)));
                        }
                        // update total factor with speceies factor
                        if (stoichiometry == 1) {
                            factorExpr = Expression.mult(factorExpr, speciesFactor);
                        } else if (stoichiometry > 1) {
                            // rxnProbExpr * (structSize^stoichiometry)
                            Expression powerExpr = Expression.power(speciesFactor, new Expression(stoichiometry));
                            factorExpr = Expression.mult(factorExpr, powerExpr);
                        }
                        if (rxnProbabilityExpr == null) {
                            rxnProbabilityExpr = new Expression(tempExpr);
                        } else {
                            rxnProbabilityExpr = Expression.mult(rxnProbabilityExpr, tempExpr);
                        }
                    }
                }
            }
        }
        // Now construct the probability expression.
        if (rateConstantExpr == null) {
            throw new MappingException("Can not find reaction rate constant in reaction: " + reactionStep.getName());
        } else if (rxnProbabilityExpr == null) {
            probExp = new Expression(rateConstantExpr);
        } else if ((rateConstantExpr != null) && (rxnProbabilityExpr != null)) {
            probExp = Expression.mult(rateConstantExpr, rxnProbabilityExpr);
        }
        // simplify the factor
        RationalExp factorRatExp = RationalExpUtils.getRationalExp(factorExpr);
        factorExpr = new Expression(factorRatExp.infixString());
        factorExpr.bindExpression(this);
        // get probability rate with converting factor
        probExp = Expression.mult(probExp, factorExpr);
        probExp = probExp.flatten();
    // //
    // // round trip to rational expression for simplifying terms like KMOLE/KMOLE ...
    // // we don't want to loose the symbol binding ... so we make a temporary symbolTable from the original binding.
    // //
    // final Expression finalExp = new Expression(probExp);
    // SymbolTable symbolTable = new SymbolTable(){
    // public void getEntries(Map<String, SymbolTableEntry> entryMap) {
    // throw new RuntimeException("should not be called");
    // }
    // public SymbolTableEntry getEntry(String identifierString) throws ExpressionBindingException {
    // return finalExp.getSymbolBinding(identifierString);
    // }
    // };
    // cbit.vcell.matrix.RationalExp ratExp = cbit.vcell.parser.RationalExpUtils.getRationalExp(probExp);
    // probExp = new Expression(ratExp.infixString());
    // probExp.bindExpression(symbolTable);
    } catch (ExpressionException e) {
        e.printStackTrace();
    }
    return probExp;
}

Example 60 with Membrane

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

the class StructureAnalyzer method refreshTotalMatrices.

/**
 * This method was created in VisualAge.
 */
private void refreshTotalMatrices() throws Exception {
    // System.out.println("StructureAnalyzer.refreshTotalMatrices()");
    // 
    // update scheme matrix for full system (slow and fast)
    // 
    ReactionSpec[] reactionSpecs = new ReactionSpec[reactionSteps.length];
    for (int j = 0; j < reactionSteps.length; j++) {
        reactionSpecs[j] = mathMapping_4_8.getSimulationContext().getReactionContext().getReactionSpec(reactionSteps[j]);
    }
    // 
    // initialize rate expressions for speciesContext's due to scheme matrix
    // 
    totalSchemeMatrix = new RationalNumberMatrix(speciesContextMappings.length, reactionSteps.length);
    for (int i = 0; i < speciesContextMappings.length; i++) {
        SpeciesContextMapping scm = speciesContextMappings[i];
        SpeciesContext sc = scm.getSpeciesContext();
        // 
        // collect slow rate expression (fast handled by FastSystem)
        // 
        Expression exp = new Expression(0.0);
        for (int j = 0; j < reactionSteps.length; j++) {
            int stoichiometry = reactionSteps[j].getStoichiometry(sc);
            totalSchemeMatrix.set_elem(i, j, stoichiometry);
            if (stoichiometry != 0) {
                if (!(reactionSteps[j] instanceof DummyReactionStep) && !reactionSpecs[j].isFast() && !reactionSpecs[j].isExcluded()) {
                    ReactionParticipant[] rps1 = reactionSteps[j].getReactionParticipants();
                    ReactionParticipant rp0 = null;
                    for (ReactionParticipant rp : rps1) {
                        if (rp.getSpeciesContext() == sc) {
                            rp0 = rp;
                            break;
                        }
                    }
                    Structure structure = reactionSteps[j].getStructure();
                    // 
                    if (rp0 != null) {
                        Expression reactRateExp = getReactionRateExpression(reactionSteps[j], rp0).renameBoundSymbols(mathMapping_4_8.getNameScope());
                        if ((structure instanceof Membrane) && (sc.getStructure() != structure)) {
                            Membrane membrane = (Membrane) structure;
                            MembraneMapping membraneMapping = (MembraneMapping) mathMapping_4_8.getSimulationContext().getGeometryContext().getStructureMapping(membrane);
                            Parameter fluxCorrectionParameter = mathMapping_4_8.getFluxCorrectionParameter(membraneMapping, (Feature) sc.getStructure());
                            Expression fluxCorrection = new Expression(fluxCorrectionParameter, mathMapping_4_8.getNameScope());
                            if (reactionSteps[j] instanceof FluxReaction) {
                                exp = Expression.add(exp, Expression.mult(fluxCorrection, reactRateExp));
                            // Expression.add(exp,new Expression(fluxCorrectionParameterSymbolName+"*"+expInfix));
                            } else if (reactionSteps[j] instanceof SimpleReaction) {
                                ModelUnitSystem unitSystem = mathMapping_4_8.getSimulationContext().getModel().getUnitSystem();
                                Expression unitFactor = mathMapping_4_8.getUnitFactor(unitSystem.getVolumeSubstanceUnit().divideBy(unitSystem.getMembraneSubstanceUnit()));
                                exp = Expression.add(exp, Expression.mult(fluxCorrection, unitFactor, reactRateExp));
                            // exp = Expression.add(exp,new Expression(fluxCorrectionParameterSymbolName+"*"+ReservedSymbol.KMOLE.getName()+"*"+expInfix));
                            } else {
                                throw new RuntimeException("Internal Error: expected ReactionStep " + reactionSteps[j] + " to be of type SimpleReaction or FluxReaction");
                            }
                        } else {
                            exp = Expression.add(exp, reactRateExp);
                        }
                    }
                }
            }
        }
        // exp.bindExpression(mathMapping);
        scm.setRate(exp.flatten());
    }
    // 
    if (totalSchemeMatrix.getNumRows() > 1) {
        totalNullSpaceMatrix = (RationalMatrix) totalSchemeMatrix.findNullSpace();
    } else {
        totalNullSpaceMatrix = null;
    }
// if (totalNullSpaceMatrix==null){
// System.out.println("total system has full rank");
// }else{
// System.out.println("StructureAnalyzer.refreshTotalMatrices(), nullSpace matrix:");
// totalNullSpaceMatrix.show();
// }
}