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

Example 6 with MassActionKinetics

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

the class TransformMassActionTableModel method saveTransformedReactions.

public void saveTransformedReactions() throws Exception {
    // isTransformable and TransformedREactions are stored according to the indexes in model reaction steps.
    boolean[] isTransformable = getTransformMassActions().getIsTransformable();
    TransformMassActions.TransformedReaction[] trs = getTransformMassActions().getTransformedReactionSteps();
    ReactionStep[] origReactions = getModel().getReactionSteps();
    // names of those who can be transformed and will be transformed
    String okTransReacNames = "";
    // names of those who can be transformed and will not be transformed
    String noTransReacNames = "";
    // names of those who can not be transformed
    String errReacNames = "";
    for (int i = 0; i < isTransformable.length; i++) {
        if (trs[i].getTransformType() == TransformMassActions.TransformedReaction.TRANSFORMABLE && getIsSelected(i)) {
            okTransReacNames = okTransReacNames + origReactions[i].getName() + ",";
        } else if (trs[i].getTransformType() == TransformMassActions.TransformedReaction.TRANSFORMABLE && !getIsSelected(i)) {
            noTransReacNames = noTransReacNames + origReactions[i].getName() + ",";
        } else if (!isTransformable[i]) {
            errReacNames = errReacNames + origReactions[i].getName() + ",";
        }
    }
    // set transformed Mass Action kinetics to model reactions
    for (int i = 0; i < origReactions.length; i++) {
        if (getIsSelected(i)) {
            // for simple reaction, we replace the original kinetics with MassActionKinetics if it wasn't MassActionKinetics
            if (origReactions[i] instanceof SimpleReaction) {
                if (!(origReactions[i].getKinetics() instanceof MassActionKinetics)) {
                    // ***Below we will physically change the simple reaction***
                    // put all kinetic parameters together into array newKps
                    Vector<Kinetics.KineticsParameter> newKps = new Vector<Kinetics.KineticsParameter>();
                    // get original kinetic parameters which are not current density and reaction rate.
                    // those parameters are basically the symbols in rate expression.
                    Vector<Kinetics.KineticsParameter> origKps = new Vector<Kinetics.KineticsParameter>();
                    Kinetics.KineticsParameter[] Kps = origReactions[i].getKinetics().getKineticsParameters();
                    for (int j = 0; j < Kps.length; j++) {
                        if (!(Kps[j].getRole() == Kinetics.ROLE_CurrentDensity || Kps[j].getRole() == Kinetics.ROLE_ReactionRate)) {
                            origKps.add(Kps[j]);
                        }
                    }
                    // create mass action kinetics for the original reaction step
                    MassActionKinetics maKinetics = new MassActionKinetics(origReactions[i]);
                    maKinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward).setExpression(trs[i].getMassActionFunction().getForwardRate());
                    maKinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse).setExpression(trs[i].getMassActionFunction().getReverseRate());
                    // Kinetics)
                    for (int j = 0; j < maKinetics.getKineticsParameters().length; j++) {
                        newKps.add(maKinetics.getKineticsParameters(j));
                    }
                    // copy other kinetic parameters from original kinetics
                    for (int j = 0; j < origKps.size(); j++) {
                        newKps.add(origKps.elementAt(j));
                    }
                    // add parameters to mass action kinetics
                    KineticsParameter[] newParameters = new KineticsParameter[newKps.size()];
                    newParameters = (KineticsParameter[]) newKps.toArray(newParameters);
                    maKinetics.addKineticsParameters(newParameters);
                    // after adding all the parameters, we bind the forward/reverse rate expression with symbol table (the reaction step itself)
                    origReactions[i].getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_KForward).getExpression().bindExpression(origReactions[i]);
                    origReactions[i].getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_KReverse).getExpression().bindExpression(origReactions[i]);
                }
            }
        // for flux, we set the flux reaction back, coz we will parse it to mass action form in stochastic math mapping.
        // However, we don't physically change it.
        }
    }
    String msg = "";
    if (!okTransReacNames.equals("")) {
        msg = msg + okTransReacNames + " have been transformed.\n";
    }
    // message to be displayed in popupdialog of DocumentWindow
    if (!msg.equals("")) {
        throw new Exception(msg);
    }
}

Example 7 with MassActionKinetics

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

the class NetworkTransformer method transform.

private void transform(SimulationContext simContext, SimulationContext transformedSimulationContext, ArrayList<ModelEntityMapping> entityMappings, MathMappingCallback mathMappingCallback, NetworkGenerationRequirements networkGenerationRequirements) {
    String msg = "Generating network: flattening...";
    mathMappingCallback.setMessage(msg);
    TaskCallbackMessage tcm = new TaskCallbackMessage(TaskCallbackStatus.Clean, "");
    simContext.appendToConsole(tcm);
    tcm = new TaskCallbackMessage(TaskCallbackStatus.TaskStart, msg);
    simContext.appendToConsole(tcm);
    long startTime = System.currentTimeMillis();
    System.out.println("Convert to bngl, execute BNG, retrieve the results.");
    try {
        BNGOutputSpec outputSpec = generateNetwork(simContext, mathMappingCallback, networkGenerationRequirements);
        if (mathMappingCallback.isInterrupted()) {
            msg = "Canceled by user.";
            tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
            simContext.appendToConsole(tcm);
            throw new UserCancelException(msg);
        }
        long endTime = System.currentTimeMillis();
        long elapsedTime = endTime - startTime;
        System.out.println("     " + elapsedTime + " milliseconds");
        Model model = transformedSimulationContext.getModel();
        ReactionContext reactionContext = transformedSimulationContext.getReactionContext();
        // ---- Parameters -----------------------------------------------------------------------------------------------
        startTime = System.currentTimeMillis();
        for (int i = 0; i < outputSpec.getBNGParams().length; i++) {
            BNGParameter p = outputSpec.getBNGParams()[i];
            // System.out.println(i+1 + ":\t\t"+ p.toString());
            if (model.getRbmModelContainer().getParameter(p.getName()) != null) {
                // if it's already there we don't try to add it again; this should be true for all of them!
                continue;
            }
            String s = p.getName();
            FakeSeedSpeciesInitialConditionsParameter fakeICParam = FakeSeedSpeciesInitialConditionsParameter.fromString(s);
            if (speciesEquivalenceMap.containsKey(fakeICParam)) {
                // we get rid of the fake parameters we use as keys
                continue;
            }
            FakeReactionRuleRateParameter fakeKineticParam = FakeReactionRuleRateParameter.fromString(s);
            if (fakeKineticParam != null) {
                System.out.println("found fakeKineticParam " + fakeKineticParam.fakeParameterName);
                // we get rid of the fake parameters we use as keys
                continue;
            }
            throw new RuntimeException("unexpected parameter " + p.getName() + " in internal BNG processing");
        // Expression exp = new Expression(p.getValue());
        // exp.bindExpression(model.getRbmModelContainer().getSymbolTable());
        // model.getRbmModelContainer().addParameter(p.getName(), exp, model.getUnitSystem().getInstance_TBD());
        }
        endTime = System.currentTimeMillis();
        elapsedTime = endTime - startTime;
        msg = "Adding " + outputSpec.getBNGParams().length + " parameters to model, " + elapsedTime + " ms";
        System.out.println(msg);
        // ---- Species ------------------------------------------------------------------------------------------------------------
        mathMappingCallback.setMessage("generating network: adding species...");
        mathMappingCallback.setProgressFraction(progressFractionQuota / 4.0f);
        startTime = System.currentTimeMillis();
        System.out.println("\nSpecies :");
        // the reactions will need this map to recover the names of species knowing only the networkFileIndex
        HashMap<Integer, String> speciesMap = new HashMap<Integer, String>();
        LinkedHashMap<String, Species> sMap = new LinkedHashMap<String, Species>();
        LinkedHashMap<String, SpeciesContext> scMap = new LinkedHashMap<String, SpeciesContext>();
        LinkedHashMap<String, BNGSpecies> crossMap = new LinkedHashMap<String, BNGSpecies>();
        List<SpeciesContext> noMapForThese = new ArrayList<SpeciesContext>();
        // final int decimalTickCount = Math.max(outputSpec.getBNGSpecies().length/10, 1);
        for (int i = 0; i < outputSpec.getBNGSpecies().length; i++) {
            BNGSpecies s = outputSpec.getBNGSpecies()[i];
            // System.out.println(i+1 + ":\t\t"+ s.toString());
            String key = s.getConcentration().infix();
            FakeSeedSpeciesInitialConditionsParameter fakeParam = FakeSeedSpeciesInitialConditionsParameter.fromString(key);
            if (fakeParam != null) {
                Pair<SpeciesContext, Expression> value = speciesEquivalenceMap.get(fakeParam);
                // the species context of the original model
                SpeciesContext originalsc = value.one;
                Expression initial = value.two;
                // replace the fake initial condition with the real one
                s.setConcentration(initial);
                // we'll have to find the species context from the cloned model which correspond to the original species
                SpeciesContext sc = model.getSpeciesContext(originalsc.getName());
                // System.out.println(sc.getName() + ", " + sc.getSpecies().getCommonName() + "   ...is one of the original seed species.");
                // existing name
                speciesMap.put(s.getNetworkFileIndex(), sc.getName());
                sMap.put(sc.getName(), sc.getSpecies());
                scMap.put(sc.getName(), sc);
                crossMap.put(sc.getName(), s);
                noMapForThese.add(sc);
                continue;
            }
            // all these species are new!
            // generate unique name for the species
            int count = 0;
            String speciesName = null;
            String nameRoot = "s";
            String speciesPatternNameString = s.extractName();
            while (true) {
                speciesName = nameRoot + count;
                if (Model.isNameUnused(speciesName, model) && !sMap.containsKey(speciesName) && !scMap.containsKey(speciesName)) {
                    break;
                }
                count++;
            }
            // newly created name
            speciesMap.put(s.getNetworkFileIndex(), speciesName);
            SpeciesContext speciesContext;
            if (s.hasCompartment()) {
                String speciesPatternCompartmentString = s.extractCompartment();
                speciesContext = new SpeciesContext(new Species(speciesName, s.getName()), model.getStructure(speciesPatternCompartmentString), null);
            } else {
                speciesContext = new SpeciesContext(new Species(speciesName, s.getName()), model.getStructure(0), null);
            }
            speciesContext.setName(speciesName);
            try {
                if (speciesPatternNameString != null) {
                    SpeciesPattern sp = RbmUtils.parseSpeciesPattern(speciesPatternNameString, model);
                    speciesContext.setSpeciesPattern(sp);
                }
            } catch (ParseException e) {
                e.printStackTrace();
                throw new RuntimeException("Bad format for species pattern string: " + e.getMessage());
            }
            // speciesContext.setSpeciesPatternString(speciesPatternString);
            // model.addSpecies(speciesContext.getSpecies());
            // model.addSpeciesContext(speciesContext);
            sMap.put(speciesName, speciesContext.getSpecies());
            scMap.put(speciesName, speciesContext);
            crossMap.put(speciesName, s);
            // }
            if (mathMappingCallback.isInterrupted()) {
                msg = "Canceled by user.";
                tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
                simContext.appendToConsole(tcm);
                throw new UserCancelException(msg);
            }
        // if(i%50 == 0) {
        // System.out.println(i+"");
        // }
        // if(i%decimalTickCount == 0) {
        // int multiplier = i/decimalTickCount;
        // float progress = progressFractionQuota/4.0f + progressFractionQuotaSpecies*multiplier;
        // mathMappingCallback.setProgressFraction(progress);
        // }
        }
        for (SpeciesContext sc1 : model.getSpeciesContexts()) {
            boolean found = false;
            for (Map.Entry<String, SpeciesContext> entry : scMap.entrySet()) {
                SpeciesContext sc2 = entry.getValue();
                if (sc1.getName().equals(sc2.getName())) {
                    found = true;
                    // System.out.println("found species context " + sc1.getName() + " of species " + sc1.getSpecies().getCommonName() + " // " + sc2.getSpecies().getCommonName());
                    break;
                }
            }
            if (found == false) {
                // we add to the map the species context and the species which exist in the model but which are not in the map yet
                // the only ones in this situation should be plain species which were not given to bngl for flattening (they are flat already)
                // System.out.println("species context " + sc1.getName() + " not found in the map. Adding it.");
                scMap.put(sc1.getName(), sc1);
                sMap.put(sc1.getName(), sc1.getSpecies());
                noMapForThese.add(sc1);
            }
        }
        for (Species s1 : model.getSpecies()) {
            boolean found = false;
            for (Map.Entry<String, Species> entry : sMap.entrySet()) {
                Species s2 = entry.getValue();
                if (s1.getCommonName().equals(s2.getCommonName())) {
                    found = true;
                    // System.out.println("found species " + s1.getCommonName());
                    break;
                }
            }
            if (found == false) {
                System.err.println("species " + s1.getCommonName() + " not found in the map!");
            }
        }
        SpeciesContext[] sca = new SpeciesContext[scMap.size()];
        scMap.values().toArray(sca);
        Species[] sa = new HashSet<Species>(sMap.values()).toArray(new Species[0]);
        model.setSpecies(sa);
        model.setSpeciesContexts(sca);
        boolean isSpatial = transformedSimulationContext.getGeometry().getDimension() > 0;
        for (SpeciesContext sc : sca) {
            if (noMapForThese.contains(sc)) {
                continue;
            }
            SpeciesContextSpec scs = reactionContext.getSpeciesContextSpec(sc);
            Parameter param = scs.getParameter(SpeciesContextSpec.ROLE_InitialConcentration);
            BNGSpecies s = crossMap.get(sc.getName());
            param.setExpression(s.getConcentration());
            SpeciesContext origSpeciesContext = simContext.getModel().getSpeciesContext(s.getName());
            if (origSpeciesContext != null) {
                ModelEntityMapping em = new ModelEntityMapping(origSpeciesContext, sc);
                entityMappings.add(em);
            } else {
                ModelEntityMapping em = new ModelEntityMapping(new GeneratedSpeciesSymbolTableEntry(sc), sc);
                if (isSpatial) {
                    scs.initializeForSpatial();
                }
                entityMappings.add(em);
            }
        }
        // for(SpeciesContext sc : sca) {		// clean all the species patterns from the flattened species, we have no sp now
        // sc.setSpeciesPattern(null);
        // }
        endTime = System.currentTimeMillis();
        elapsedTime = endTime - startTime;
        msg = "Adding " + outputSpec.getBNGSpecies().length + " species to model, " + elapsedTime + " ms";
        System.out.println(msg);
        // ---- Reactions -----------------------------------------------------------------------------------------------------
        mathMappingCallback.setMessage("generating network: adding reactions...");
        mathMappingCallback.setProgressFraction(progressFractionQuota / 4.0f * 3.0f);
        startTime = System.currentTimeMillis();
        System.out.println("\nReactions :");
        Map<String, HashSet<String>> ruleKeyMap = new HashMap<String, HashSet<String>>();
        Map<String, BNGReaction> directBNGReactionsMap = new HashMap<String, BNGReaction>();
        Map<String, BNGReaction> reverseBNGReactionsMap = new HashMap<String, BNGReaction>();
        for (int i = 0; i < outputSpec.getBNGReactions().length; i++) {
            BNGReaction r = outputSpec.getBNGReactions()[i];
            if (!r.isRuleReversed()) {
                // direct
                directBNGReactionsMap.put(r.getKey(), r);
            } else {
                reverseBNGReactionsMap.put(r.getKey(), r);
            }
            // 
            // for each rule name, store set of keySets (number of unique keysets are number of generated reactions from this ruleName).
            // 
            HashSet<String> keySet = ruleKeyMap.get(r.getRuleName());
            if (keySet == null) {
                keySet = new HashSet<String>();
                ruleKeyMap.put(r.getRuleName(), keySet);
            }
            keySet.add(r.getKey());
        }
        Map<String, ReactionStep> reactionStepMap = new HashMap<String, ReactionStep>();
        for (int i = 0; i < outputSpec.getBNGReactions().length; i++) {
            BNGReaction bngReaction = outputSpec.getBNGReactions()[i];
            // System.out.println(i+1 + ":\t\t"+ r.writeReaction());
            String baseName = bngReaction.getRuleName();
            String reactionName = null;
            HashSet<String> keySetsForThisRule = ruleKeyMap.get(bngReaction.getRuleName());
            if (keySetsForThisRule.size() == 1 && model.getReactionStep(bngReaction.getRuleName()) == null && !reactionStepMap.containsKey(bngReaction.getRuleName())) {
                // we can reuse the reaction rule labels
                reactionName = bngReaction.getRuleName();
            } else {
                reactionName = bngReaction.getRuleName() + "_0";
                while (true) {
                    if (model.getReactionStep(reactionName) == null && !reactionStepMap.containsKey(reactionName)) {
                        // we can reuse the reaction rule labels
                        break;
                    }
                    reactionName = TokenMangler.getNextEnumeratedToken(reactionName);
                }
            }
            // 
            if (directBNGReactionsMap.containsValue(bngReaction)) {
                BNGReaction forwardBNGReaction = bngReaction;
                BNGReaction reverseBNGReaction = reverseBNGReactionsMap.get(bngReaction.getKey());
                String name = forwardBNGReaction.getRuleName();
                if (name.endsWith(ReactionRule.DirectHalf)) {
                    name = name.substring(0, name.indexOf(ReactionRule.DirectHalf));
                }
                if (name.endsWith(ReactionRule.InverseHalf)) {
                    name = name.substring(0, name.indexOf(ReactionRule.InverseHalf));
                }
                ReactionRule rr = model.getRbmModelContainer().getReactionRule(name);
                Structure structure = rr.getStructure();
                boolean bReversible = reverseBNGReaction != null;
                SimpleReaction sr = new SimpleReaction(model, structure, reactionName, bReversible);
                for (int j = 0; j < forwardBNGReaction.getReactants().length; j++) {
                    BNGSpecies s = forwardBNGReaction.getReactants()[j];
                    String scName = speciesMap.get(s.getNetworkFileIndex());
                    SpeciesContext sc = model.getSpeciesContext(scName);
                    Reactant reactant = sr.getReactant(scName);
                    if (reactant == null) {
                        int stoichiometry = 1;
                        sr.addReactant(sc, stoichiometry);
                    } else {
                        int stoichiometry = reactant.getStoichiometry();
                        stoichiometry += 1;
                        reactant.setStoichiometry(stoichiometry);
                    }
                }
                for (int j = 0; j < forwardBNGReaction.getProducts().length; j++) {
                    BNGSpecies s = forwardBNGReaction.getProducts()[j];
                    String scName = speciesMap.get(s.getNetworkFileIndex());
                    SpeciesContext sc = model.getSpeciesContext(scName);
                    Product product = sr.getProduct(scName);
                    if (product == null) {
                        int stoichiometry = 1;
                        sr.addProduct(sc, stoichiometry);
                    } else {
                        int stoichiometry = product.getStoichiometry();
                        stoichiometry += 1;
                        product.setStoichiometry(stoichiometry);
                    }
                }
                MassActionKinetics targetKinetics = new MassActionKinetics(sr);
                sr.setKinetics(targetKinetics);
                KineticsParameter kforward = targetKinetics.getForwardRateParameter();
                KineticsParameter kreverse = targetKinetics.getReverseRateParameter();
                String kforwardNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionForwardRate).getName();
                if (!kforward.getName().equals(kforwardNewName)) {
                    targetKinetics.renameParameter(kforward.getName(), kforwardNewName);
                    kforward = targetKinetics.getForwardRateParameter();
                }
                final String kreverseNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionReverseRate).getName();
                if (!kreverse.getName().equals(kreverseNewName)) {
                    targetKinetics.renameParameter(kreverse.getName(), kreverseNewName);
                    kreverse = targetKinetics.getReverseRateParameter();
                }
                applyKineticsExpressions(forwardBNGReaction, kforward, targetKinetics);
                if (reverseBNGReaction != null) {
                    applyKineticsExpressions(reverseBNGReaction, kreverse, targetKinetics);
                }
                // String fieldParameterName = kforward.getName();
                // fieldParameterName += "_" + r.getRuleName();
                // kforward.setName(fieldParameterName);
                reactionStepMap.put(reactionName, sr);
            } else if (reverseBNGReactionsMap.containsValue(bngReaction) && !directBNGReactionsMap.containsKey(bngReaction.getKey())) {
                // reverse only (must be irreversible)
                BNGReaction reverseBNGReaction = reverseBNGReactionsMap.get(bngReaction.getKey());
                ReactionRule rr = model.getRbmModelContainer().getReactionRule(reverseBNGReaction.extractRuleName());
                Structure structure = rr.getStructure();
                boolean bReversible = false;
                SimpleReaction sr = new SimpleReaction(model, structure, reactionName, bReversible);
                for (int j = 0; j < reverseBNGReaction.getReactants().length; j++) {
                    BNGSpecies s = reverseBNGReaction.getReactants()[j];
                    String scName = speciesMap.get(s.getNetworkFileIndex());
                    SpeciesContext sc = model.getSpeciesContext(scName);
                    Reactant reactant = sr.getReactant(scName);
                    if (reactant == null) {
                        int stoichiometry = 1;
                        sr.addReactant(sc, stoichiometry);
                    } else {
                        int stoichiometry = reactant.getStoichiometry();
                        stoichiometry += 1;
                        reactant.setStoichiometry(stoichiometry);
                    }
                }
                for (int j = 0; j < reverseBNGReaction.getProducts().length; j++) {
                    BNGSpecies s = reverseBNGReaction.getProducts()[j];
                    String scName = speciesMap.get(s.getNetworkFileIndex());
                    SpeciesContext sc = model.getSpeciesContext(scName);
                    Product product = sr.getProduct(scName);
                    if (product == null) {
                        int stoichiometry = 1;
                        sr.addProduct(sc, stoichiometry);
                    } else {
                        int stoichiometry = product.getStoichiometry();
                        stoichiometry += 1;
                        product.setStoichiometry(stoichiometry);
                    }
                }
                MassActionKinetics k = new MassActionKinetics(sr);
                sr.setKinetics(k);
                KineticsParameter kforward = k.getForwardRateParameter();
                KineticsParameter kreverse = k.getReverseRateParameter();
                String kforwardNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionForwardRate).getName();
                if (!kforward.getName().equals(kforwardNewName)) {
                    k.renameParameter(kforward.getName(), kforwardNewName);
                    kforward = k.getForwardRateParameter();
                }
                final String kreverseNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionReverseRate).getName();
                if (!kreverse.getName().equals(kreverseNewName)) {
                    k.renameParameter(kreverse.getName(), kreverseNewName);
                    kreverse = k.getReverseRateParameter();
                }
                applyKineticsExpressions(reverseBNGReaction, kforward, k);
                // String fieldParameterName = kforward.getName();
                // fieldParameterName += "_" + r.getRuleName();
                // kforward.setName(fieldParameterName);
                reactionStepMap.put(reactionName, sr);
            }
        }
        for (ReactionStep rs : model.getReactionSteps()) {
            reactionStepMap.put(rs.getName(), rs);
        }
        ReactionStep[] reactionSteps = new ReactionStep[reactionStepMap.size()];
        reactionStepMap.values().toArray(reactionSteps);
        model.setReactionSteps(reactionSteps);
        if (mathMappingCallback.isInterrupted()) {
            msg = "Canceled by user.";
            tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
            simContext.appendToConsole(tcm);
            throw new UserCancelException(msg);
        }
        endTime = System.currentTimeMillis();
        elapsedTime = endTime - startTime;
        msg = "Adding " + outputSpec.getBNGReactions().length + " reactions to model, " + elapsedTime + " ms";
        System.out.println(msg);
        // clean all the reaction rules
        model.getRbmModelContainer().getReactionRuleList().clear();
        // ---- Observables -------------------------------------------------------------------------------------------------
        mathMappingCallback.setMessage("generating network: adding observables...");
        mathMappingCallback.setProgressFraction(progressFractionQuota / 8.0f * 7.0f);
        startTime = System.currentTimeMillis();
        System.out.println("\nObservables :");
        RbmModelContainer rbmmc = model.getRbmModelContainer();
        for (int i = 0; i < outputSpec.getObservableGroups().length; i++) {
            ObservableGroup o = outputSpec.getObservableGroups()[i];
            if (rbmmc.getParameter(o.getObservableGroupName()) != null) {
                System.out.println("   ...already exists.");
                // if it's already there we don't try to add it again; this should be true for all of them!
                continue;
            }
            ArrayList<Expression> terms = new ArrayList<Expression>();
            for (int j = 0; j < o.getListofSpecies().length; j++) {
                Expression term = Expression.mult(new Expression(o.getSpeciesMultiplicity()[j]), new Expression(speciesMap.get(o.getListofSpecies()[j].getNetworkFileIndex())));
                terms.add(term);
            }
            Expression exp = Expression.add(terms.toArray(new Expression[terms.size()])).flatten();
            exp.bindExpression(rbmmc.getSymbolTable());
            RbmObservable originalObservable = rbmmc.getObservable(o.getObservableGroupName());
            VCUnitDefinition observableUnitDefinition = originalObservable.getUnitDefinition();
            rbmmc.removeObservable(originalObservable);
            Parameter newParameter = rbmmc.addParameter(o.getObservableGroupName(), exp, observableUnitDefinition);
            RbmObservable origObservable = simContext.getModel().getRbmModelContainer().getObservable(o.getObservableGroupName());
            ModelEntityMapping em = new ModelEntityMapping(origObservable, newParameter);
            entityMappings.add(em);
        }
        if (mathMappingCallback.isInterrupted()) {
            msg = "Canceled by user.";
            tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
            simContext.appendToConsole(tcm);
            throw new UserCancelException(msg);
        }
        endTime = System.currentTimeMillis();
        elapsedTime = endTime - startTime;
        msg = "Adding " + outputSpec.getObservableGroups().length + " observables to model, " + elapsedTime + " ms";
        System.out.println(msg);
    } catch (PropertyVetoException ex) {
        ex.printStackTrace(System.out);
        throw new RuntimeException(ex.getMessage());
    } catch (ExpressionBindingException ex) {
        ex.printStackTrace(System.out);
        throw new RuntimeException(ex.getMessage());
    } catch (ModelException ex) {
        ex.printStackTrace(System.out);
        throw new RuntimeException(ex.getMessage());
    } catch (ExpressionException ex) {
        ex.printStackTrace(System.out);
        throw new RuntimeException(ex.getMessage());
    } catch (ClassNotFoundException ex) {
        throw new RuntimeException(ex.getMessage());
    } catch (IOException ex) {
        throw new RuntimeException(ex.getMessage());
    }
    System.out.println("Done transforming");
    msg = "Generating math...";
    System.out.println(msg);
    mathMappingCallback.setMessage(msg);
    mathMappingCallback.setProgressFraction(progressFractionQuota);
}

Example 8 with MassActionKinetics

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

the class RulebasedTransformer method transform.

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

Example 9 with MassActionKinetics

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

the class ITextWriter method writeReactions.

// each reaction has its own table, ordered by the structures.
protected void writeReactions(Chapter physioChapter, Model model) throws DocumentException {
    if (model == null) {
        return;
    }
    Paragraph reactionParagraph = new Paragraph();
    reactionParagraph.add(new Chunk("Structures and Reactions Diagram").setLocalDestination(model.getName()));
    Section reactionDiagramSection = physioChapter.addSection(reactionParagraph, physioChapter.numberDepth() + 1);
    try {
        addImage(reactionDiagramSection, encodeJPEG(generateDocReactionsImage(model, null)));
    } catch (Exception e) {
        e.printStackTrace();
        throw new DocumentException(e.getClass().getName() + ": " + e.getMessage());
    }
    for (int i = 0; i < model.getNumStructures(); i++) {
        ReactionStep[] reactionSteps = model.getReactionSteps();
        ReactionStep rs = null;
        Table modifierTable = null;
        Table reactionTable = null;
        boolean firstTime = true;
        Section reactStructSection = null;
        for (int j = 0; j < reactionSteps.length; j++) {
            if (reactionSteps[j].getStructure() == model.getStructure(i)) {
                // can also use structureName1.equals(structureName2)
                if (firstTime) {
                    Paragraph linkParagraph = new Paragraph();
                    linkParagraph.add(new Chunk("Reaction(s) in " + model.getStructure(i).getName()).setLocalDestination(model.getStructure(i).getName()));
                    reactStructSection = physioChapter.addSection(linkParagraph, physioChapter.numberDepth() + 1);
                    firstTime = false;
                }
                rs = reactionSteps[j];
                String type;
                if (rs instanceof SimpleReaction) {
                    type = "Reaction";
                } else {
                    type = "Flux";
                }
                // write Reaction equation as a table
                // Get the image arrow cell depending on type of reactionStep : MassAction => double arrow, otherwise, forward arrow
                boolean bReversible = false;
                if (rs.getKinetics() instanceof MassActionKinetics) {
                    bReversible = true;
                }
                Cell arrowImageCell = getReactionArrowImageCell(bReversible);
                // Get reactants and products strings
                ReactionCanvas rc = new ReactionCanvas();
                rc.setReactionStep(rs);
                ReactionCanvasDisplaySpec rcdSpec = rc.getReactionCanvasDisplaySpec();
                String reactants = rcdSpec.getLeftText();
                String products = rcdSpec.getRightText();
                // Create table and add cells for reactants, arrow(s) images, products
                int[] widths = { 8, 1, 8 };
                reactionTable = getTable(3, 100, 0, 2, 2);
                // Add reactants as cell
                Cell tableCell = createCell(reactants, getBold());
                tableCell.setHorizontalAlignment(Cell.ALIGN_RIGHT);
                tableCell.setBorderColor(Color.white);
                reactionTable.addCell(tableCell);
                // add arrow(s) image as cell
                if (arrowImageCell != null) {
                    arrowImageCell.setHorizontalAlignment(Cell.ALIGN_CENTER);
                    arrowImageCell.setBorderColor(Color.white);
                    reactionTable.addCell(arrowImageCell);
                }
                // add products as cell
                tableCell = createCell(products, getBold());
                tableCell.setBorderColor(Color.white);
                reactionTable.addCell(tableCell);
                // reactionTable.setBorderColor(Color.white);
                reactionTable.setWidths(widths);
                // Identify modifiers,
                ReactionParticipant[] rpArr = rs.getReactionParticipants();
                Vector<ReactionParticipant> modifiersVector = new Vector<ReactionParticipant>();
                for (int k = 0; k < rpArr.length; k += 1) {
                    if (rpArr[k] instanceof Catalyst) {
                        modifiersVector.add(rpArr[k]);
                    }
                }
                // Write the modifiers in a separate table, if present
                if (modifiersVector.size() > 0) {
                    modifierTable = getTable(1, 50, 0, 1, 1);
                    modifierTable.addCell(createCell("Modifiers List", getBold(DEF_HEADER_FONT_SIZE), 1, 1, Element.ALIGN_CENTER, true));
                    StringBuffer modifierNames = new StringBuffer();
                    for (int k = 0; k < modifiersVector.size(); k++) {
                        modifierNames.append(((Catalyst) modifiersVector.elementAt(k)).getName() + "\n");
                    }
                    modifierTable.addCell(createCell(modifierNames.toString().trim(), getFont()));
                    modifiersVector.removeAllElements();
                }
                Section reactionSection = reactStructSection.addSection(type + " " + rs.getName(), reactStructSection.numberDepth() + 1);
                // Annotation
                VCMetaData vcMetaData = rs.getModel().getVcMetaData();
                if (vcMetaData.getFreeTextAnnotation(rs) != null) {
                    Table annotTable = getTable(1, 100, 1, 3, 3);
                    annotTable.addCell(createCell("Reaction Annotation", getBold(DEF_HEADER_FONT_SIZE), 1, 1, Element.ALIGN_CENTER, true));
                    annotTable.addCell(createCell(vcMetaData.getFreeTextAnnotation(rs), getFont()));
                    reactionSection.add(annotTable);
                // reactionSection.add(new Paragraph("\""+rs.getAnnotation()+"\""));
                }
                // reaction table
                if (reactionTable != null) {
                    reactionSection.add(reactionTable);
                    // re-set reactionTable
                    reactionTable = null;
                }
                if (modifierTable != null) {
                    reactionSection.add(modifierTable);
                    modifierTable = null;
                }
                // Write kinetics parameters, etc. in a table
                writeKineticsParams(reactionSection, rs);
            }
        }
    }
}

Example 10 with MassActionKinetics

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

the class ApplicationConstraintsGenerator method steadyStateFromApplication.

/**
 * Insert the method's description here.
 * Creation date: (6/26/01 8:25:55 AM)
 * @return cbit.vcell.constraints.ConstraintContainerImpl
 */
public static ConstraintContainerImpl steadyStateFromApplication(SimulationContext simContext, double tolerance) {
    try {
        ConstraintContainerImpl ccImpl = new ConstraintContainerImpl();
        // ====================
        // add physical limits
        // ====================
        // 
        // no negative concentrations
        // 
        cbit.vcell.model.Model model = simContext.getModel();
        cbit.vcell.model.SpeciesContext[] speciesContexts = model.getSpeciesContexts();
        for (int i = 0; i < speciesContexts.length; i++) {
            ccImpl.addSimpleBound(new SimpleBounds(speciesContexts[i].getName(), new RealInterval(0, Double.POSITIVE_INFINITY), AbstractConstraint.PHYSICAL_LIMIT, "non-negative concentration"));
        }
        for (int i = 0; i < speciesContexts.length; i++) {
            SpeciesContextSpecParameter initParam = (simContext.getReactionContext().getSpeciesContextSpec(speciesContexts[i])).getInitialConditionParameter();
            if (initParam != null) {
                double initialValue = initParam.getExpression().evaluateConstant();
                double lowInitialValue = Math.min(initialValue / tolerance, initialValue * tolerance);
                double highInitialValue = Math.max(initialValue / tolerance, initialValue * tolerance);
                ccImpl.addSimpleBound(new SimpleBounds(speciesContexts[i].getName(), new RealInterval(lowInitialValue, highInitialValue), AbstractConstraint.MODELING_ASSUMPTION, "close to specified \"initialCondition\""));
            }
        }
        // =========================
        // add modeling assumptions
        // =========================
        // 
        // mass action forward and reverse rates should be non-negative
        // 
        cbit.vcell.model.ReactionStep[] reactionSteps = model.getReactionSteps();
        for (int i = 0; i < reactionSteps.length; i++) {
            Kinetics kinetics = reactionSteps[i].getKinetics();
            if (kinetics instanceof MassActionKinetics) {
                Expression forwardRateConstraintExp = new Expression(((MassActionKinetics) kinetics).getForwardRateParameter().getExpression().infix() + ">=0");
                forwardRateConstraintExp = getSteadyStateExpression(forwardRateConstraintExp);
                if (!forwardRateConstraintExp.compareEqual(new Expression(1.0))) {
                    ccImpl.addGeneralConstraint(new GeneralConstraint(forwardRateConstraintExp, AbstractConstraint.MODELING_ASSUMPTION, "non-negative forward rate"));
                }
                Expression reverseRateConstraintExp = new Expression(((MassActionKinetics) kinetics).getReverseRateParameter().getExpression().infix() + ">=0");
                reverseRateConstraintExp = getSteadyStateExpression(reverseRateConstraintExp);
                if (!reverseRateConstraintExp.compareEqual(new Expression(1.0))) {
                    ccImpl.addGeneralConstraint(new GeneralConstraint(reverseRateConstraintExp, AbstractConstraint.MODELING_ASSUMPTION, "non-negative reverse rate"));
                }
            }
            KineticsParameter authoritativeParameter = kinetics.getAuthoritativeParameter();
            Expression kineticRateConstraintExp = new Expression(authoritativeParameter.getName() + "==" + authoritativeParameter.getExpression().infix());
            kineticRateConstraintExp = getSteadyStateExpression(kineticRateConstraintExp);
            if (!kineticRateConstraintExp.compareEqual(new Expression(1.0))) {
                ccImpl.addGeneralConstraint(new GeneralConstraint(kineticRateConstraintExp, AbstractConstraint.MODELING_ASSUMPTION, "definition"));
            }
        }
        // 
        try {
            simContext.setMathDescription(simContext.createNewMathMapping().getMathDescription());
        } catch (Throwable e) {
            e.printStackTrace(System.out);
            throw new RuntimeException("cannot create mathDescription");
        }
        MathDescription mathDesc = simContext.getMathDescription();
        if (mathDesc.getGeometry().getDimension() > 0) {
            throw new RuntimeException("spatial simulations not yet supported");
        }
        CompartmentSubDomain subDomain = (CompartmentSubDomain) mathDesc.getSubDomains().nextElement();
        java.util.Enumeration<Equation> enumEquations = subDomain.getEquations();
        while (enumEquations.hasMoreElements()) {
            Equation equation = (Equation) enumEquations.nextElement();
            Expression rateConstraintExp = new Expression(equation.getRateExpression().infix() + "==0");
            rateConstraintExp = getSteadyStateExpression(rateConstraintExp);
            if (!rateConstraintExp.compareEqual(new Expression(1.0))) {
                // not a trivial constraint (always true)
                ccImpl.addGeneralConstraint(new GeneralConstraint(rateConstraintExp, AbstractConstraint.PHYSICAL_LIMIT, "definition of steady state"));
            }
        }
        // 
        for (int i = 0; i < reactionSteps.length; i++) {
            Kinetics kinetics = reactionSteps[i].getKinetics();
            Kinetics.KineticsParameter[] parameters = kinetics.getKineticsParameters();
            for (int j = 0; j < parameters.length; j++) {
                Expression exp = parameters[j].getExpression();
                if (exp.getSymbols() == null || exp.getSymbols().length == 0) {
                    // 
                    try {
                        double constantValue = exp.evaluateConstant();
                        double lowValue = Math.min(constantValue / tolerance, constantValue * tolerance);
                        double highValue = Math.max(constantValue / tolerance, constantValue * tolerance);
                        RealInterval interval = new RealInterval(lowValue, highValue);
                        ccImpl.addSimpleBound(new SimpleBounds(parameters[j].getName(), interval, AbstractConstraint.MODELING_ASSUMPTION, "parameter close to model default"));
                    } catch (cbit.vcell.parser.ExpressionException e) {
                        System.out.println("error evaluating parameter " + parameters[j].getName() + " in reaction step " + reactionSteps[i].getName());
                    }
                } else {
                    Expression parameterDefinitionExp = new Expression(parameters[j].getName() + "==" + parameters[j].getExpression().infix());
                    ccImpl.addGeneralConstraint(new GeneralConstraint(getSteadyStateExpression(parameterDefinitionExp), AbstractConstraint.MODELING_ASSUMPTION, "parameter definition"));
                }
            }
        }
        ccImpl.addSimpleBound(new SimpleBounds(model.getFARADAY_CONSTANT().getName(), new RealInterval(model.getFARADAY_CONSTANT().getExpression().evaluateConstant()), AbstractConstraint.PHYSICAL_LIMIT, "Faraday's constant"));
        ccImpl.addSimpleBound(new SimpleBounds(model.getTEMPERATURE().getName(), new RealInterval(300), AbstractConstraint.PHYSICAL_LIMIT, "Absolute Temperature Kelvin"));
        ccImpl.addSimpleBound(new SimpleBounds(model.getGAS_CONSTANT().getName(), new RealInterval(model.getGAS_CONSTANT().getExpression().evaluateConstant()), AbstractConstraint.PHYSICAL_LIMIT, "ideal gas constant"));
        ccImpl.addSimpleBound(new SimpleBounds(model.getKMILLIVOLTS().getName(), new RealInterval(model.getKMILLIVOLTS().getExpression().evaluateConstant()), AbstractConstraint.PHYSICAL_LIMIT, "ideal gas constant"));
        // 
        // add K_fluxs
        // 
        java.util.Enumeration<Variable> enumVars = mathDesc.getVariables();
        while (enumVars.hasMoreElements()) {
            Variable var = (Variable) enumVars.nextElement();
            if (var.getName().startsWith("Kflux_") && var instanceof Function) {
                Expression kfluxExp = new Expression(((Function) var).getExpression());
                kfluxExp.bindExpression(mathDesc);
                kfluxExp = MathUtilities.substituteFunctions(kfluxExp, mathDesc);
                kfluxExp = kfluxExp.flatten();
                ccImpl.addSimpleBound(new SimpleBounds(var.getName(), new RealInterval(kfluxExp.evaluateConstant()), AbstractConstraint.MODELING_ASSUMPTION, "flux conversion factor"));
            }
        }
        return ccImpl;
    } catch (cbit.vcell.parser.ExpressionException e) {
        e.printStackTrace(System.out);
        return null;
    } catch (java.beans.PropertyVetoException e) {
        e.printStackTrace(System.out);
        return null;
    }
}