use of cbit.vcell.model.Feature in project vcell by virtualcell.
the class FRAPStudy method createNewRefBioModel.
public static BioModel createNewRefBioModel(FRAPStudy sourceFrapStudy, String baseDiffusionRate, TimeStep tStep, KeyValue simKey, User owner, FieldDataIdentifierSpec psfFDIS, int startingIndexForRecovery) throws Exception {
if (owner == null) {
throw new Exception("Owner is not defined");
}
ROI cellROI_2D = sourceFrapStudy.getFrapData().getRoi(FRAPData.VFRAP_ROI_ENUM.ROI_CELL.name());
Extent extent = sourceFrapStudy.getFrapData().getImageDataset().getExtent();
TimeBounds timeBounds = FRAPOptData.getEstimatedRefTimeBound(sourceFrapStudy);
double timeStepVal = FRAPOptData.REFERENCE_DIFF_DELTAT;
int numX = cellROI_2D.getRoiImages()[0].getNumX();
int numY = cellROI_2D.getRoiImages()[0].getNumY();
int numZ = cellROI_2D.getRoiImages().length;
short[] shortPixels = cellROI_2D.getRoiImages()[0].getPixels();
byte[] bytePixels = new byte[numX * numY * numZ];
final byte EXTRACELLULAR_PIXVAL = 0;
final byte CYTOSOL_PIXVAL = 1;
for (int i = 0; i < bytePixels.length; i++) {
if (shortPixels[i] != 0) {
bytePixels[i] = CYTOSOL_PIXVAL;
}
}
VCImage maskImage;
try {
maskImage = new VCImageUncompressed(null, bytePixels, extent, numX, numY, numZ);
} catch (ImageException e) {
e.printStackTrace();
throw new RuntimeException("failed to create mask image for geometry");
}
Geometry geometry = new Geometry("geometry", maskImage);
if (geometry.getGeometrySpec().getNumSubVolumes() != 2) {
throw new Exception("Cell ROI has no ExtraCellular.");
}
int subVolume0PixVal = ((ImageSubVolume) geometry.getGeometrySpec().getSubVolume(0)).getPixelValue();
geometry.getGeometrySpec().getSubVolume(0).setName((subVolume0PixVal == EXTRACELLULAR_PIXVAL ? EXTRACELLULAR_NAME : CYTOSOL_NAME));
int subVolume1PixVal = ((ImageSubVolume) geometry.getGeometrySpec().getSubVolume(1)).getPixelValue();
geometry.getGeometrySpec().getSubVolume(1).setName((subVolume1PixVal == CYTOSOL_PIXVAL ? CYTOSOL_NAME : EXTRACELLULAR_NAME));
geometry.getGeometrySurfaceDescription().updateAll();
BioModel bioModel = new BioModel(null);
bioModel.setName("unnamed");
Model model = new Model("model");
bioModel.setModel(model);
Feature extracellular = model.addFeature(EXTRACELLULAR_NAME);
Feature cytosol = model.addFeature(CYTOSOL_NAME);
Membrane plasmaMembrane = model.addMembrane(PLASMAMEMBRANE_NAME);
String roiDataName = FRAPStudy.ROI_EXTDATA_NAME;
final int ONE_DIFFUSION_SPECIES_COUNT = 1;
final int MOBILE_SPECIES_INDEX = 0;
Expression[] diffusionConstants = new Expression[ONE_DIFFUSION_SPECIES_COUNT];
Species[] species = new Species[ONE_DIFFUSION_SPECIES_COUNT];
SpeciesContext[] speciesContexts = new SpeciesContext[ONE_DIFFUSION_SPECIES_COUNT];
Expression[] initialConditions = new Expression[ONE_DIFFUSION_SPECIES_COUNT];
// Mobile Species
diffusionConstants[MOBILE_SPECIES_INDEX] = new Expression(baseDiffusionRate);
species[MOBILE_SPECIES_INDEX] = new Species(SPECIES_NAME_PREFIX_MOBILE, "Mobile bleachable species");
speciesContexts[MOBILE_SPECIES_INDEX] = new SpeciesContext(null, species[MOBILE_SPECIES_INDEX].getCommonName(), species[MOBILE_SPECIES_INDEX], cytosol);
FieldFunctionArguments postBleach_first = new FieldFunctionArguments(roiDataName, "postbleach_first", new Expression(0), VariableType.VOLUME);
FieldFunctionArguments prebleach_avg = new FieldFunctionArguments(roiDataName, "prebleach_avg", new Expression(0), VariableType.VOLUME);
Expression expPostBleach_first = new Expression(postBleach_first.infix());
Expression expPreBleach_avg = new Expression(prebleach_avg.infix());
initialConditions[MOBILE_SPECIES_INDEX] = Expression.div(expPostBleach_first, expPreBleach_avg);
SimulationContext simContext = new SimulationContext(bioModel.getModel(), geometry);
bioModel.addSimulationContext(simContext);
FeatureMapping cytosolFeatureMapping = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(cytosol);
FeatureMapping extracellularFeatureMapping = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(extracellular);
MembraneMapping plasmaMembraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(plasmaMembrane);
SubVolume cytSubVolume = geometry.getGeometrySpec().getSubVolume(CYTOSOL_NAME);
SubVolume exSubVolume = geometry.getGeometrySpec().getSubVolume(EXTRACELLULAR_NAME);
SurfaceClass pmSurfaceClass = geometry.getGeometrySurfaceDescription().getSurfaceClass(exSubVolume, cytSubVolume);
cytosolFeatureMapping.setGeometryClass(cytSubVolume);
extracellularFeatureMapping.setGeometryClass(exSubVolume);
plasmaMembraneMapping.setGeometryClass(pmSurfaceClass);
cytosolFeatureMapping.getUnitSizeParameter().setExpression(new Expression(1.0));
extracellularFeatureMapping.getUnitSizeParameter().setExpression(new Expression(1.0));
plasmaMembraneMapping.getUnitSizeParameter().setExpression(new Expression(1.0));
for (int i = 0; i < initialConditions.length; i++) {
model.addSpecies(species[i]);
model.addSpeciesContext(speciesContexts[i]);
}
for (int i = 0; i < speciesContexts.length; i++) {
SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(speciesContexts[i]);
scs.getInitialConditionParameter().setExpression(initialConditions[i]);
scs.getDiffusionParameter().setExpression(diffusionConstants[i]);
}
MathMapping mathMapping = simContext.createNewMathMapping();
MathDescription mathDesc = mathMapping.getMathDescription();
// Add PSF function
mathDesc.addVariable(new Function(Simulation.PSF_FUNCTION_NAME, new Expression(psfFDIS.getFieldFuncArgs().infix()), null));
simContext.setMathDescription(mathDesc);
SimulationVersion simVersion = new SimulationVersion(simKey, "sim1", owner, new GroupAccessNone(), new KeyValue("0"), new BigDecimal(0), new Date(), VersionFlag.Current, "", null);
Simulation newSimulation = new Simulation(simVersion, simContext.getMathDescription());
newSimulation.getSolverTaskDescription().setSolverDescription(SolverDescription.FiniteVolumeStandalone);
simContext.addSimulation(newSimulation);
newSimulation.getSolverTaskDescription().setTimeBounds(timeBounds);
newSimulation.getSolverTaskDescription().setOutputTimeSpec(new UniformOutputTimeSpec(timeStepVal));
newSimulation.getMeshSpecification().setSamplingSize(cellROI_2D.getISize());
newSimulation.getSolverTaskDescription().setTimeStep(new TimeStep(timeStepVal, timeStepVal, timeStepVal));
return bioModel;
}
use of cbit.vcell.model.Feature in project vcell by virtualcell.
the class DiffEquMathMapping method computeBoundaryConditionSource.
// determine membrane inside and outside subvolume
public static Pair<SubVolume, SubVolume> computeBoundaryConditionSource(Model model, SimulationContext simContext, SurfaceClass surfaceClass) {
SubVolume outerSubVolume = null;
SubVolume innerSubVolume = null;
Structure[] mappedStructures = simContext.getGeometryContext().getStructuresFromGeometryClass(surfaceClass);
// inside and outside correspond to structure hierarchy when present
for (Structure s : mappedStructures) {
if (s instanceof Membrane) {
Membrane m = (Membrane) s;
Feature infeature = model.getStructureTopology().getInsideFeature(m);
if (infeature != null) {
FeatureMapping insm = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(infeature);
if (insm.getGeometryClass() instanceof SubVolume) {
innerSubVolume = (SubVolume) insm.getGeometryClass();
}
}
Feature outfeature = model.getStructureTopology().getOutsideFeature(m);
if (outfeature != null) {
FeatureMapping outsm = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(outfeature);
if (outsm.getGeometryClass() instanceof SubVolume) {
outerSubVolume = (SubVolume) outsm.getGeometryClass();
}
}
}
}
// make the choice deterministic
if (innerSubVolume == null || outerSubVolume == null || innerSubVolume == outerSubVolume) {
Set<SubVolume> sv = surfaceClass.getAdjacentSubvolumes();
Iterator<SubVolume> iterator = sv.iterator();
innerSubVolume = iterator.next();
outerSubVolume = iterator.next();
if (innerSubVolume.getName().compareTo(outerSubVolume.getName()) > 0) {
SubVolume temp = innerSubVolume;
innerSubVolume = outerSubVolume;
outerSubVolume = temp;
}
}
Pair<SubVolume, SubVolume> ret = new Pair<>(innerSubVolume, outerSubVolume);
return ret;
}
use of cbit.vcell.model.Feature in project vcell by virtualcell.
the class DiffEquMathMapping 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);
}
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();
Domain domain = null;
if (sm.getGeometryClass() != null) {
domain = new Domain(sm.getGeometryClass());
}
if (struct instanceof Feature || struct instanceof Membrane) {
if (sm.getGeometryClass() instanceof SurfaceClass) {
if (scs.isWellMixed()) {
scm.setVariable(new MembraneRegionVariable(scm.getSpeciesContext().getName(), domain));
} else if (getSimulationContext().isStoch() && getSimulationContext().getGeometry().getDimension() > 0 && !scs.isForceContinuous()) {
scm.setVariable(new MemVariable(scm.getSpeciesContext().getName() + "_conc", domain));
} else {
scm.setVariable(new MemVariable(scm.getSpeciesContext().getName(), domain));
}
} else {
if (scs.isWellMixed()) {
scm.setVariable(new VolumeRegionVariable(scm.getSpeciesContext().getName(), domain));
} else if (getSimulationContext().isStoch() && getSimulationContext().getGeometry().getDimension() > 0 && !scs.isForceContinuous()) {
scm.setVariable(new VolVariable(scm.getSpeciesContext().getName() + "_conc", domain));
} else {
scm.setVariable(new VolVariable(scm.getSpeciesContext().getName(), domain));
}
}
} else {
throw new MappingException("class " + scm.getSpeciesContext().getStructure().getClass() + " not supported");
}
mathSymbolMapping.put(scm.getSpeciesContext(), scm.getVariable().getName());
}
}
}
use of cbit.vcell.model.Feature in project vcell by virtualcell.
the class StructureAnalyzer method refreshFastMatrices.
/**
* This method was created in VisualAge.
*/
private void refreshFastMatrices() throws Exception {
// System.out.println("StructureAnalyzer.refreshFastMatrices()");
//
// update scheme matrix for fast system
//
fastSchemeMatrix = new RationalNumberMatrix(fastSpeciesContextMappings.length, fastReactionSteps.length);
for (int i = 0; i < fastSpeciesContextMappings.length; i++) {
for (int j = 0; j < fastReactionSteps.length; j++) {
fastSchemeMatrix.set_elem(i, j, fastReactionSteps[j].getStoichiometry(fastSpeciesContextMappings[i].getSpeciesContext()));
}
}
//
for (int i = 0; i < fastSpeciesContextMappings.length; i++) {
SpeciesContextMapping scm = fastSpeciesContextMappings[i];
SpeciesContext sc = scm.getSpeciesContext();
//
// collect fast rate expression
//
Expression exp = new Expression(0.0);
for (int j = 0; j < fastReactionSteps.length; j++) {
int stoichiometry = fastReactionSteps[j].getStoichiometry(sc);
ReactionSpec reactionSpec = mathMapping_4_8.getSimulationContext().getReactionContext().getReactionSpec(fastReactionSteps[j]);
if (stoichiometry != 0) {
if (!reactionSpec.isFast()) {
throw new Exception("expected only fast rates");
}
if (reactionSpec.isExcluded()) {
throw new Exception("expected only included rates");
}
ReactionParticipant[] rps = fastReactionSteps[j].getReactionParticipants();
ReactionParticipant rp0 = null;
for (ReactionParticipant rp : rps) {
if (rp.getSpeciesContext() == sc) {
rp0 = rp;
break;
}
}
//
if (rp0 != null) {
Structure structure = fastReactionSteps[j].getStructure();
Expression fastRateExpression = getReactionRateExpression(fastReactionSteps[j], rp0).renameBoundSymbols(mathMapping_4_8.getNameScope());
if ((structure instanceof Membrane) && (rp0.getStructure() != structure)) {
Membrane membrane = (Membrane) structure;
MembraneMapping membraneMapping = (MembraneMapping) mathMapping_4_8.getSimulationContext().getGeometryContext().getStructureMapping(membrane);
Expression fluxCorrection = new Expression(mathMapping_4_8.getFluxCorrectionParameter(membraneMapping, (Feature) rp0.getStructure()), mathMapping_4_8.getNameScope());
exp = Expression.add(exp, Expression.mult(fluxCorrection, fastRateExpression));
} else {
exp = Expression.add(exp, new Expression(fastRateExpression));
}
}
}
}
// exp.bindExpression(mathMapping);
scm.setFastRate(exp.flatten());
}
// System.out.println("StructureAnalyzer.refreshFastMatrices(), scheme matrix:");
// fastSchemeMatrix.show();
//
// update null space matrix
//
fastNullSpaceMatrix = fastSchemeMatrix.findNullSpace();
// if (fastNullSpaceMatrix==null){
// System.out.println("fast system has full rank");
// }else{
// System.out.println("StructureAnalyzer.refreshFastMatrices(), nullSpace matrix:");
// fastNullSpaceMatrix.show();
// }
}
use of cbit.vcell.model.Feature in project vcell by virtualcell.
the class StructureAnalyzer method getReactionRateExpression.
public Expression getReactionRateExpression(ReactionStep reactionStep, ReactionParticipant reactionParticipant) throws Exception {
if (reactionParticipant instanceof Catalyst) {
throw new Exception("Catalyst " + reactionParticipant + " doesn't have a rate for this reaction");
// return new Expression(0.0);
}
double stoich = reactionStep.getStoichiometry(reactionParticipant.getSpeciesContext());
if (stoich == 0.0) {
return new Expression(0.0);
}
if (reactionStep.getKinetics() instanceof DistributedKinetics) {
DistributedKinetics distributedKinetics = (DistributedKinetics) reactionStep.getKinetics();
if (stoich != 1) {
Expression exp = Expression.mult(new Expression(stoich), new Expression(distributedKinetics.getReactionRateParameter(), mathMapping_4_8.getNameScope()));
return exp;
} else {
Expression exp = new Expression(distributedKinetics.getReactionRateParameter(), mathMapping_4_8.getNameScope());
return exp;
}
} else if (reactionStep.getKinetics() instanceof LumpedKinetics) {
Structure.StructureSize structureSize = reactionStep.getStructure().getStructureSize();
//
// need to put this into concentration/time with respect to structure for reaction.
//
LumpedKinetics lumpedKinetics = (LumpedKinetics) reactionStep.getKinetics();
Expression factor = null;
ModelUnitSystem unitSystem = mathMapping_4_8.getSimulationContext().getModel().getUnitSystem();
if (reactionStep.getStructure() instanceof Feature || ((reactionStep.getStructure() instanceof Membrane) && reactionStep instanceof FluxReaction)) {
VCUnitDefinition lumpedToVolumeSubstance = unitSystem.getVolumeSubstanceUnit().divideBy(unitSystem.getLumpedReactionSubstanceUnit());
factor = Expression.div(new Expression(lumpedToVolumeSubstance.getDimensionlessScale()), new Expression(structureSize, mathMapping_4_8.getNameScope()));
} else if (reactionStep.getStructure() instanceof Membrane && reactionStep instanceof SimpleReaction) {
VCUnitDefinition lumpedToVolumeSubstance = unitSystem.getMembraneSubstanceUnit().divideBy(unitSystem.getLumpedReactionSubstanceUnit());
factor = Expression.div(new Expression(lumpedToVolumeSubstance.getDimensionlessScale()), new Expression(structureSize, mathMapping_4_8.getNameScope()));
} else {
throw new RuntimeException("failed to create reaction rate expression for reaction " + reactionStep.getName() + ", with kinetic type of " + reactionStep.getKinetics().getClass().getName());
}
if (stoich != 1) {
Expression exp = Expression.mult(new Expression(stoich), Expression.mult(new Expression(lumpedKinetics.getLumpedReactionRateParameter(), mathMapping_4_8.getNameScope()), factor));
return exp;
} else {
Expression exp = Expression.mult(new Expression(lumpedKinetics.getLumpedReactionRateParameter(), mathMapping_4_8.getNameScope()), factor);
return exp;
}
} else {
throw new RuntimeException("unexpected kinetic type " + reactionStep.getKinetics().getClass().getName());
}
}
Aggregations