Implement chemotherapy process in v3

src/arcade/patch/agent/cell/PatchCell.java

@@ -19,6 +19,7 @@
 import arcade.patch.agent.module.PatchModuleApoptosis;
 import arcade.patch.agent.module.PatchModuleMigration;
 import arcade.patch.agent.module.PatchModuleProliferation;
+import arcade.patch.agent.process.PatchProcessChemotherapy;
 import arcade.patch.agent.process.PatchProcessMetabolism;
 import arcade.patch.agent.process.PatchProcessSignaling;
 import arcade.patch.env.grid.PatchGrid;
@@ -341,6 +342,8 @@ public Process makeProcess(ProcessDomain domain, String version) {
                 return PatchProcessMetabolism.make(this, version);
             case SIGNALING:
                 return PatchProcessSignaling.make(this, version);
+            case CHEMOTHERAPY:
+                return PatchProcessChemotherapy.make(this, version);
             case UNDEFINED:
             default:
                 return null;
@@ -396,6 +399,10 @@ public void step(SimState simstate) {
             }
         }
 
+        if (processes.get(Domain.CHEMOTHERAPY) != null) {
+            processes.get(Domain.CHEMOTHERAPY).step(simstate.random, sim);
+        }
+
         // Step the module for the cell state.
         if (module != null) {
             module.step(simstate.random, sim);

src/arcade/patch/agent/process/PatchProcessChemotherapy.java

@@ -0,0 +1,128 @@
+package arcade.patch.agent.process;
+
+import sim.util.Bag;
+import ec.util.MersenneTwisterFast;
+import arcade.core.sim.Simulation;
+import arcade.core.util.MiniBox;
+import arcade.patch.agent.cell.PatchCell;
+import arcade.patch.env.grid.PatchGrid;
+
+/**
+ * Implementation of {@link Process} for cell chemotherapy.
+ *
+ * <p>The {@code PatchProcessChemotherapy} process:
+ *
+ * <ul>
+ *   <li>gets available drugs from the environment
+ *   <li>calculates drug uptake given cell size and state
+ *   <li>steps the chemotherapy process to determine changes to internal drug concentration
+ *   <li>updates drug environment with uptake and decay
+ *   <li>determines whether cell should apoptose due to chemotherapy
+ * </ul>
+ */
+public abstract class PatchProcessChemotherapy extends PatchProcess {
+    /** Threshold at which cells undergo apoptosis. */
+    protected final double chemotherapyThreshold;
+
+    /** Constant drug uptake rate [fmol drug/um^2 cell/min/M drug]. */
+    protected final double drugUptakeRate;
+
+    /** Rate at which drugs are removed from the cell. */
+    protected final double drugRemovalRate;
+
+    /** Rate at which drugs decay in the cell. */
+    protected final double drugDecayRate;
+
+    /** Volume of cell [um<sup>3</sup>]. */
+    double volume;
+
+    /** Volume fraction. */
+    protected double f;
+
+    /** Internal amount of the drug. */
+    protected double intAmt;
+
+    /** External amout of the drug. */
+    protected double extAmt;
+
+    /** Uptake of the drug in the current step. */
+    protected double uptakeAmt;
+
+    /** Whether the cell was killed by chemotherapy. */
+    protected boolean wasChemo;
+
+    /**
+     * Creates a chemotherapy {@code Process} for the given {@link PatchCell}.
+     *
+     * @param cell the {@link PatchCell} the process is associated with
+     */
+    PatchProcessChemotherapy(PatchCell cell) {
+        super(cell);
+
+        // Initialize process
+        volume = cell.getVolume();
+
+        // Set loaded parameters.
+        MiniBox parameters = cell.getParameters();
+        chemotherapyThreshold = parameters.getDouble("chemotherapy/CHEMOTHERAPY_THRESHOLD");
+        drugUptakeRate = parameters.getDouble("chemotherapy/CONSTANT_DRUG_UPTAKE_RATE");
+        drugRemovalRate = parameters.getDouble("chemotherapy/DRUG_REMOVAL_RATE");
+        drugDecayRate = parameters.getDouble("chemotherapy/DRUG_DECAY_RATE");
+
+        // Initial internal concentrations.
+        extAmt = 0.0;
+        uptakeAmt = 0.0;
+        intAmt = 0.0;
+    }
+
+    /**
+     * Steps the chemotherapy process.
+     *
+     * @param random the random number generator
+     * @param sim the simulation instance
+     */
+    abstract void stepProcess(MersenneTwisterFast random, Simulation sim);
+
+    /**
+     * Gets the external drug concentrations from the environment.
+     *
+     * @param sim the simulation instance
+     */
+    private void updateExternal(Simulation sim) {
+        extAmt = sim.getLattice("DRUG").getAverageValue(location) * location.getVolume();
+        extAmt *= (1.0 - drugDecayRate);
+    }
+
+    @Override
+    public void step(MersenneTwisterFast random, Simulation sim) {
+        // Calculate fraction of volume occupied by cell.
+        Bag bag = ((PatchGrid) sim.getGrid()).getObjectsAtLocation(location);
+        double totalVolume = PatchCell.calculateTotalVolume(bag);
+        f = volume / totalVolume;
+
+        // Get external drug concentration.
+        updateExternal(sim);
+
+        // Calculate drug uptake and internal concentration.
+        stepProcess(random, sim);
+
+        // Update environment for the drug.
+        sim.getLattice("DRUG").updateValue(location, 1.0 - uptakeAmt / extAmt);
+    }
+
+    /**
+     * Creates a {@code PatchProcessChemotherapy} for the given version.
+     *
+     * @param cell the {@link PatchCell} the process is associated with
+     * @param version the process version
+     * @return the process instance
+     */
+    public static PatchProcess make(PatchCell cell, String version) {
+        switch (version.toUpperCase()) {
+            case "SIMPLE":
+                return new PatchProcessChemotherapySimple(cell);
+            default:
+                return null;
+        }
+    }
+}

src/arcade/patch/agent/process/PatchProcessChemotherapySimple.java

@@ -0,0 +1,75 @@
+package arcade.patch.agent.process;
+
+import ec.util.MersenneTwisterFast;
+import arcade.core.agent.process.Process;
+import arcade.core.sim.Simulation;
+import arcade.patch.agent.cell.PatchCell;
+import static arcade.patch.util.PatchEnums.State;
+
+/**
+ * Extension of {@link PatchProcessChemotherapy} for simple chemotherapy.
+ *
+ * <p>{@code PatchProcessChemotherapySimple} assumes a constant drug uptake rate and a threshold for
+ * apoptosis.
+ */
+public class PatchProcessChemotherapySimple extends PatchProcessChemotherapy {
+    /**
+     * Creates a simple chemotherapy {@code Process} for the given {@link PatchCell}.
+     *
+     * <p>Loaded parameters include:
+     *
+     * <ul>
+     *   <li>{@code CONSTANT_DRUG_UPTAKE_RATE} = constant drug uptake rate
+     *   <li>{@code KILL_THRESHOLD} = drug kill threshold concentration
+     *   <li>{@code KILL_RATE} = constant kill rate
+     * </ul>
+     *
+     * @param cell the {@link PatchCell} the process is associated with
+     */
+    public PatchProcessChemotherapySimple(PatchCell cell) {
+        super(cell);
+    }
+
+    @Override
+    public void stepProcess(MersenneTwisterFast random, Simulation sim) {
+        double drugInt = intAmt;
+        double drugExt = extAmt;
+
+        // Calculate drug uptake rate based on concentration gradient.
+        double area = location.getArea() * f;
+        double surfaceArea = area * 2 + (volume / area) * location.getPerimeter(f);
+        double drugGrad = (extAmt / location.getVolume()) - (drugInt / volume);
+        drugGrad *= drugGrad < 1E-10 ? 0 : 1;
+        double drugUptake = drugUptakeRate * drugGrad * surfaceArea;
+        drugInt += drugUptake;
+
+        // If drug concentration exceeds kill threshold kill cells with probability based on drug
+        // concentration.
+        if (cell.getState() == State.PROLIFERATIVE && drugInt > chemotherapyThreshold) {
+            double oxygen = sim.getLattice("OXYGEN").getAverageValue(location);
+            double p = Math.pow(oxygen, 2) / (Math.pow(oxygen, 2) + Math.pow(drugInt, 2));
+
+            // TODO: Update probability
+            if (random.nextDouble() < p) {
+                cell.setState(State.APOPTOTIC);
+                wasChemo = true;
+            }
+        }
+        intAmt = Math.exp(-drugRemovalRate) * drugInt;
+        // intAmt = drugInt;
+        uptakeAmt = drugUptake;
+    }
+
+    @Override
+    public void update(Process process) {
+        PatchProcessChemotherapySimple chemotherapy = (PatchProcessChemotherapySimple) process;
+        double split = this.cell.getVolume() / this.volume;
+
+        // Update this process as split of given process.
+        this.volume = this.cell.getVolume();
+        this.intAmt = chemotherapy.intAmt * split;
+
+        chemotherapy.volume = chemotherapy.cell.getVolume();
+        chemotherapy.intAmt *= (1 - split);
+    }
+}

src/arcade/patch/env/component/PatchComponentDegrade.java

@@ -75,7 +75,6 @@ public void schedule(Schedule schedule) {
     @Override
     public void register(Simulation sim, String layer) {
         Component component = sim.getComponent(layer);
-
         if (!(component instanceof PatchComponentSitesGraph)) {
             return;
         }

src/arcade/patch/env/component/PatchComponentDose.java

@@ -0,0 +1,122 @@
+package arcade.patch.env.component;
+
+import java.util.ArrayList;
+import sim.engine.Schedule;
+import sim.engine.SimState;
+import arcade.core.env.component.Component;
+import arcade.core.env.lattice.Lattice;
+import arcade.core.env.operation.Operation;
+import arcade.core.sim.Series;
+import arcade.core.sim.Simulation;
+import arcade.core.util.MiniBox;
+import arcade.patch.env.operation.PatchOperationGenerator;
+import arcade.patch.sim.PatchSeries;
+import static arcade.patch.env.component.PatchComponentSites.SiteLayer;
+import static arcade.patch.util.PatchEnums.Category;
+import static arcade.patch.util.PatchEnums.Ordering;
+
+public class PatchComponentDose implements Component {
+    private final ArrayList<DoseLayer> layers;
+    private final int latticeHeight;
+    private final int latticeLength;
+    private final int latticeWidth;
+    private final double mediaAmount;
+    private final double mediaVolume;
+    private final double latticePatchVolume;
+    private final double latticePatchArea;
+    private final double doseStart;
+    private final double doseDuration;
+    private final double doseInterval;
+    private final double doseEnd;
+
+    public PatchComponentDose(Series series, MiniBox parameters) {
+        layers = new ArrayList<>();
+
+        latticeLength = series.length;
+        latticeWidth = series.width;
+        latticeHeight = series.height;
+
+        // Set loaded parameters.
+        mediaAmount = parameters.getDouble("MEDIA_AMOUNT");
+        doseStart = parameters.getDouble("DOSE_START");
+        doseDuration = parameters.getDouble("DOSE_DURATION");
+        doseInterval = parameters.getDouble("DOSE_INTERVAL");
+        doseEnd = parameters.getDouble("DOSE_END");
+
+        // Set patch parameters.
+        MiniBox patch = ((PatchSeries) series).patch;
+        latticePatchVolume = patch.getDouble("LATTICE_VOLUME");
+        latticePatchArea = patch.getDouble("LATTICE_AREA");
+
+        mediaVolume = latticePatchArea * latticeLength * latticeWidth * mediaAmount;
+    }
+
+    protected static class DoseLayer {
+        final String name;
+        final double[][][] current;
+        final SiteLayer siteLayer;
+        final double initialConcentration;
+        double currentAmount;
+
+        DoseLayer(String name, SiteLayer siteLayer, PatchOperationGenerator generator) {
+            this.name = name;
+            this.siteLayer = siteLayer;
+            current = generator.latticeCurrent;
+            initialConcentration = generator.concentration;
+            currentAmount = 0;
+        }
+    }
+
+    @Override
+    public void schedule(Schedule schedule) {
+        schedule.scheduleRepeating(doseStart, Ordering.FIRST_COMPONENT.ordinal(), this);
+    }
+
+    @Override
+    public void register(Simulation sim, String layer) {
+        String[] layerSplit = layer.split(":");
+        Lattice lattice = sim.getLattice(layerSplit[1]);
+        Operation generator = lattice.getOperation(Category.GENERATOR);
+        Component component = sim.getComponent(layerSplit[0]);
+
+        if (!(component instanceof PatchComponentSitesSource)) {
+            return;
+        }
+
+        PatchComponentSitesSource sites = (PatchComponentSitesSource) component;
+        SiteLayer siteLayer =
+                sites.layers.stream()
+                        .filter(sl -> sl.name.equalsIgnoreCase(layerSplit[1]))
+                        .findFirst()
+                        .orElse(null);
+
+        if (siteLayer != null) {
+            DoseLayer doseLayer =
+                    new DoseLayer(layer, siteLayer, (PatchOperationGenerator) generator);
+            layers.add(doseLayer);
+        }
+    }
+
+    @Override
+    public void step(SimState simstate) {
+        double tick = simstate.schedule.getTime();
+
+        double doseAmount = 400000;
+        for (DoseLayer layer : layers) {
+            if (tick > doseEnd) {
+                layer.currentAmount = 0;
+                layer.siteLayer.concentration = 0;
+                return;
+            }
+
+            if (tick >= doseStart && (tick - doseStart) % doseInterval < doseDuration) {
+            layer.currentAmount += doseAmount;
+            } else {
+                layer.currentAmount = 0;
+            }
+
+            layer.siteLayer.concentration = layer.currentAmount / mediaVolume;
+
+        }
+    }
+}

src/arcade/patch/env/component/PatchComponentPulse.java

@@ -137,7 +137,6 @@ public void register(Simulation sim, String layer) {
         Lattice lattice = sim.getLattice(layerSplit[1]);
         Operation generator = lattice.getOperation(Category.GENERATOR);
         Component component = sim.getComponent(layerSplit[0]);
-
         if (!(component instanceof PatchComponentSitesSource)) {
             return;
         }
@@ -174,10 +173,15 @@ public void step(SimState simstate) {
                 }
             }
 
+            System.out.println(delta);
+            System.out.println(layer.currentAmount);
             // Update available concentrations.
             layer.currentAmount = Math.max(0, layer.currentAmount - delta);
-            layer.siteLayer.concentration = layer.currentAmount / mediaVolume;
+            System.out.println(layer.currentAmount);
 
+            System.out.println(layer.siteLayer.concentration);
+            layer.siteLayer.concentration = layer.currentAmount / mediaVolume;
+            System.out.println(layer.siteLayer.concentration);
             // Pulse returns concentration to initial value.
             if (tick % pulseInterval == 0) {
                 layer.currentAmount = layer.initialConcentration * mediaVolume;