AI Screening for Immunogenic Cell Death Inducers

Analysis by the aitrendblend editorial team · Published research explained ·

• immunogenic cell death
• high throughput screening
• DIC microscopy
• model assisted labeling
• transfer learning
• YOLOv9
• DAMPs
• cold tumors

Why cold tumors need a better screening method

Immunotherapy has improved outcomes in melanoma, lymphoma and lung cancer, yet it works unevenly. The authors single out cold tumors as the stubborn case. These tumors have low immunogenicity and limited immune cell activity, so checkpoint drugs have little to switch back on. Zhang and colleagues described the goal of turning cold tumors hot in a 2022 review in Trends in Immunology, and Chen and Mellman laid out the wider logic of the cancer immunity set point in Nature in 2017. Both are cited by the paper and both frame the same ambition. Give the immune system a reason to notice the tumor.

One way to supply that reason is immunogenic cell death, or ICD. A cell dying this way does not slip away quietly. It swells, its membrane ruptures, and it releases damage associated molecular patterns, usually shortened to DAMPs. Three of them dominate this study. Calreticulin moves to the cell surface, ATP leaks into the surrounding fluid, and the chromatin protein HMGB1 escapes from the nucleus. Together they invite dendritic cells and macrophages to engulf the dying cell and present its antigens to T cells.

Finding drugs that cause this is where the practical trouble starts. Existing screens, as the authors describe them, usually pick one DAMP and read it with a colorimetric or luminescence assay. Articles in Methods in Enzymology by Zhao, Liu, Kepp and Kroemer on HMGB1 release, and by Dubyak on luciferase based ATP detection, show how those assays work. Each one captures a single molecule at a single moment. ICD is a sequence, as Galluzzi and colleagues stress in a 2023 review in the Journal of Translational Medicine. Signaling, membrane changes, rupture and release unfold over hours, and one endpoint number cannot see that sequence.

Imaging can see it, but imaging creates its own bottleneck, because someone has to look at thousands of frames. This paper sets out to remove that person from the loop, at least for the first pass. If you follow the medical imaging coverage on this site, the pattern will feel familiar. A model replaces a tedious human reading task. What is unusual is that the images are not of patients. They show cells in a dish, and the end product is a ranked list of drug candidates, not a diagnosis. You can find how we group related work in our hub for AI in medical imaging and diagnostic research.

How the pipeline works

The system has three stations. A live cell imager records differential interference contrast images, known as DIC, of cells exposed to a chemical library. A detector reads those images and counts two kinds of cell. A validation stage then tests the flagged compounds with conventional biology. Only the middle station is machine learning, and it is where the interesting decisions sit.

Borrowing labels from fluorescence

DIC is a label free technique. It shows cell shape and edges without any dye, which makes it attractive for screening and also hard to annotate, because swollen cells and healthy cells look alike to an untrained eye. The authors solved the labeling problem with two fluorescent reporters that exist only while training data is built. They engineered MDA-MB-231 breast cancer cells to carry HMGB1 tagged with turbo GFP, and in a healthy cell the green signal sits in the nucleus. They also added SYTOX Orange to the medium, a DNA dye that cannot cross an intact membrane.

That setup tells a clean three stage story. In a living cell the green stays in the nucleus and no red appears. In what the authors call a swelling cell, or SC, the cell has puffed up and the green signal has spread through the cytoplasm, but the membrane still holds, so the red dye stays out. In a dead cell, or DC, the membrane has failed, the green protein has escaped, and the red dye lights up the nucleus. Fluorescence therefore defines the ground truth, and the model is later asked to recognize the same states from DIC alone.

Model assisted labeling

Annotating swollen cells by hand is slow, and the paper puts numbers on that. Three experts needed between 4.32 and 4.66 seconds per cell. The authors built a loop instead. A small hand labeled set trains a YOLOv9 detector, the detector proposes labels on new images, experts correct the proposals, and the corrected labels return to training. Then the cycle repeats. In their comparison this model assisted labeling, or MAL, took 1.04 seconds per cell on ten microscopy images spanning the DIC, RFP and GFP channels. The model proposed 1560 labels, experts removed 70, and 1490 cells remained.

Anyone who has assembled a nuclei detection dataset will see why this pattern appeals. Reviewing a proposal is faster than drawing from scratch, and the model’s mistakes steer experts toward the hard cases.

Table 1 of the paper, reformatted. Expert scores use the MAL labels as ground truth.

The comparison with manual labeling deserves a careful read, though. Expert recall lands between 43.2 and 65.5 percent, and expert SC precision falls as low as 48.6 percent. Those figures describe how far each expert sits from the MAL set. They do not say who is right, because the MAL set has no independent check in the paper, and its own precision and recall cells in the table are blank. The experts also disagree with one another, which supports the idea that SC is a fuzzy category. Still, a sentence such as experts miss half the cells is a statement about agreement and not about accuracy. The paper also leaves open whether the human annotators could see the fluorescence channels that the MAL labels may have drawn on.

Fluorescence pretraining, then DIC only fine tuning

The second idea is what the authors call FPDF. The network is first pretrained on merged images that stack the DIC, RFP and GFP channels, so it learns what swelling and rupture look like when fluorescence makes them obvious. It is then fine tuned on DIC alone. That keeps most of the learned features while removing any need for reporters at screening time. It is a classic transfer learning move with an unusual twist, since the source and the target are the same cells seen through different channels.

Training used a large YOLOv9 variant on a single RTX 3090 with five fold cross validation, 100 epochs per fold, batch size 4 and 1024 by 1024 inputs. The optimizer was SGD at a learning rate of 0.01 with momentum 0.937 and weight decay 0.0005, and mosaic augmentation was set to 1.0 with mixup at 0.15. The paper says mosaic was switched off after the 15th epoch. Many YOLO recipes turn it off for the final 15 epochs instead, so anyone reproducing the work should check the released code to see which was meant.

A debris filter built from two numbers

Raw SC detections included debris. The authors noticed that false SC boxes were often too small or oddly proportioned, so they fed each box’s width and height into a gradient boosted tree classifier, XGBoost, to decide between SC and debris. Equation 1 gives the rule.

Here \(x_i=(w_i,h_i)\) holds the box width and height, \(f_k\) is the k th tree, \(\eta\) is the learning rate, \(\sigma\) is the sigmoid and \(\theta\) is a threshold such as 0.5. It is a tiny model with two features and a cutoff. Fig. 4B in the paper shows why it can work. Plotting width against height gives two clouds, a compact one for debris in the lower left and a broad one for SC. The classes barely overlap, and the authors report that all six classifiers they tried performed well. Those were Gaussian naive Bayes, logistic regression, SVM, random forest, LightGBM and XGBoost.

What the detection numbers actually say

The paper scores detection with precision, recall, F1 and mAP50, the mean over classes of the average precision when a predicted box must overlap the true box by at least half.

Table 2 of the paper, reformatted. SC means swelling cell and DC means dead cell.

Start with how the text describes the gains. The fine tuned model’s SC F1 of 0.470 is called about 2.6 percent better than the DIC only baseline at 0.444, and the final system’s 0.509 is called about 6.5 percent better. Those are absolute differences in score, so they are percentage points. In relative terms, which is how most readers will hear the word percent, the final SC F1 is roughly 15 percent higher, and SC mAP50 climbs from 0.255 to 0.349, roughly 37 percent. That is my arithmetic on the paper’s own table. Either way the improvement is real.

The more revealing column is precision. The DIC only baseline has the best SC precision in the table at 0.468, and the final system sits at 0.434. For DC the baseline’s 0.810 slips to 0.793. Everything the pipeline gained came from recall. SC recall rises from 0.422 to 0.615, and DC recall from 0.705 to 0.795. That is a defensible trade for a screen, since a missed swelling cell can hide a hit. But the text credits the debris filter with improving SC precision, and the table shows a gain of just 0.010, from 0.424 to 0.434, with recall unchanged at 0.615.

What does a precision of 0.434 mean in practice? Fewer than half of the SC boxes the final system draws match an annotated swelling cell at the 50 percent overlap threshold. The debris filter removed few of the wrong ones. Because it sees only width and height, it can only catch boxes of the wrong size, and Fig. 4B suggests those were already a small share. My reading is that most remaining false SC boxes are real cells caught in an in between state, the same ambiguity that made three experts disagree about them. That is a hypothesis and not a finding from the paper.

So why does a noisy detector still produce a useful screen? Because the quantity the screen reads is a curve and not a box. For each compound the system plots SC and DC counts across about 144 frames. If false SC detections hover at a roughly constant level, they raise the baseline without creating a peak. In Fig. 5A the vehicle controls, which are naive, DPBS and DMSO, look flat near zero for both classes. A real swelling wave then stands out as a transient bump. The logic is plausible and the results are consistent with it, yet the paper never tests it directly, for example by adding synthetic noise or by relating peak calls to per cell accuracy.

There is also a question about the evaluation protocol. Table 2 does not say whether the scores come from held out folds, and it reports no spread across the five folds. The ensemble pools models that were trained on different splits. If any of those scores were computed on data a member had already seen, the ensemble rows would be flattered. The released code on GitHub should settle that, and it is worth opening before you trust the ensemble numbers.

The blind test

The test used eight compounds at three concentrations of 1, 10 and 100 μM for 48 hours, with a frame every 20 minutes. The set mixed drugs reported to induce ICD with drugs without known ICD activity. The authors read the AI generated curves in Fig. 5 and saw three patterns. Flat curves meant no effect at that dose. A rise in DC with no real SC bump meant cells died without a visible swelling phase. A rise in DC accompanied by a clear SC peak meant the swelling and rupture sequence they associate with ICD.

The calls came from those peaks. SC peaks appeared for compound B at 10 and 100 μM, compound D at 100 μM, compound G at 100 μM and compound H at 1, 10 and 100 μM. The paper gives no numeric rule for what counts as an SC peak. Section 3.3 reads as though the curves were inspected, and the interpretation algorithm in Fig. 1B is drawn as a box and not specified.

Compiled from Table 3 and Fig. 6A of the paper.

Read the table from the right and the AI flagged four of eight compounds. Three of them, erastin, FIN56 and RSL3, held up as ICD inducers through every follow up assay. The fourth, raptinal, showed cleaved caspase 3 and PARP1 on the western blot, so it also triggers apoptosis, although Vernon and colleagues reported in 2022 in Molecular Cancer Research that it can drive gasdermin E dependent pyroptosis. The abstract says the AI identified three ICD inducers from eight candidates. Section 3.3 says it classified four, and the authors then treat raptinal as a mixed case that shows the AI’s sensitivity. Both readings fit the data, but the headline count depends on how you score raptinal. The four compounds the AI did not flag all behaved as plain apoptosis inducers, so the negatives look sound.

Now look at what the three clean hits share. Erastin, FIN56 and RSL3 are all ferroptosis inducers, a form of iron dependent cell death. And RSL3, compound H, is the compound the authors name in their limitations section as the source of the training data, alongside the MDA-MB-231 line. Whatever the blinding protocol was, the model had already seen what RSL3 does to these cells. The paper does not say whether the training wells and the blind test wells came from separate experiments, which would help, yet the compound identity still overlaps.

That reframes what the test proves. It shows the detector can tell a swell then rupture pattern from a plain death pattern, and that it can do so for two further ferroptosis inducers plus an unexpected pyroptosis case. It does not show that the detector can find ICD driven by a different mechanism. The authors acknowledge the narrow training data themselves. In my view the blind test framing deserves that caveat up front and not in the discussion.

capture key cellular patterns that may even elude expert analysisKim et al., Computers in Biology and Medicine, 2025, describing the system’s promise

Readers who follow ICD research may also ask about compounds A and E. Doxorubicin and oxaliplatin appear in the paper as plain apoptosis inducers with no ICD activity. Elsewhere in the field, whether drugs of those classes can provoke ICD has been discussed, and this paper does not engage with the question [RELATED STUDY NEEDED — owner to supply]. The AI and the western blots agree that A and E are not swelling type killers under these conditions, but the single cell line and the 48 hour window may matter.

What the validation assays confirm

After the AI calls, the team ran conventional assays. A crystal violet viability test gave an IC50 for each compound, from 0.0581 μM for RSL3 to 37.33 μM for etoposide, and those values set the doses for later work. Next came the apoptosis exclusion step. Live caspase 3 and 7 imaging flagged compounds C and F, while the western blot flagged A, C, E, F and G. The two assays disagree on A, E and G, and the paper moves ahead with the western blot result.

For the four candidates, LDH release and silver stained supernatant proteins both pointed to membrane rupture. Calreticulin exposure on the cell surface rose for all four. HMGB1 fell in the cell lysate and rose in the medium for all four. Extracellular ATP is where the story thins. The text says eATP levels in the supernatants were elevated. The figure shows that only RSL3 reached significance, at the three star level, while erastin, FIN56 and raptinal are marked not significant. The measurements were also taken at different times, 1, 4, 2 and 6 hours for B, D, G and H, after pretreatment with ARL 67156.

Phagocytosis was tested with mouse bone marrow derived dendritic cells and macrophages fed human cancer cells treated with each compound. Reading the bars, fold changes run roughly two to four for dendritic cells and five to seven for macrophages. The representative flow plots put absolute uptake at 1.31 percent for control dendritic cells and 0.36 percent for control macrophages, rising to between about 3.4 and 5.5 percent after treatment. Those are real signals, yet small in absolute terms.

Compound B’s dose also changes between experiments. It is 5 μM for LDH, calreticulin and phagocytosis, 50 μM for silver staining and eATP, and 10 μM in the HMGB1 legend, against an IC50 of 2.764 μM. Because the authors chose IC50 based doses to equalize cellular response, such shifts matter. Section 2.2.6 and the Fig. 7B legend also disagree about the silver staining dose, 50 versus 5 μM. It is probably a typo, but it should be fixed. The statistics are unpaired Student t tests on three replicates against a control, and the paper does not mention a correction for multiple comparisons.

None of that overturns the main result. Membrane rupture, calreticulin exposure, HMGB1 release and phagocytosis all point the same way for B, D, G and H. It does mean that the claim that all four release the full DAMP set rests mostly on calreticulin and HMGB1, while eATP supports it for RSL3 alone.

From a dish to a clinic

Clinical translation gap

This is a drug discovery tool and not a patient facing one, so the translation gap is about candidates and not about diagnoses. The distance is still long. The work uses one human breast cancer cell line grown in a plate. The immune readout relies on mouse bone marrow cells, which means human tumor cells are fed to mouse phagocytes. That pairing shows that dying cells are eaten more often, but it says little about how a human immune system would respond. Fig. 1C draws antigen presentation and T cell activation as the endpoint, yet the paper reports no T cell data and no animal experiment that tests whether treated cells protect against tumor growth. A compound that looks like an ICD inducer in this pipeline is a lead worth further study. It is not evidence of a therapy.

Regulatory and safety notes

Nothing in the paper makes a clinical or regulatory claim. The animal work, which supplied the bone marrow cells, was approved by the institutional animal care committee at Inha University. Several screened compounds, such as doxorubicin and puromycin, are cytotoxic agents used here as laboratory reagents, and the paper reports no safety or dosing data in living organisms. Treat every result as preclinical.

Limitations, with the real numbers

Sample size is the first constraint. The labeling comparison used ten images and 1490 cells. The blind test used eight compounds. The validation assays used three replicates each. The paper does not report the size of the final training set in the main text, and it reports no standard deviation across the five folds, so readers cannot tell how stable the detector is.

Dataset bias is the second. The training images came from one cell line and one ferroptosis inducer, and the authors say so in their discussion. A headline number can also hide a weak class, as we saw in our analysis of a skin lesion classifier. Here the equivalent is the SC class, with an F1 of 0.509 beside a DC F1 of 0.794.

Generalization is the third. Morphology varies by cell type, by plate coating, by imaging system and by how long cells sit in medium. The model has not been tested on another cell line, another microscope or another family of compounds. Tiny effective samples behind large augmented datasets have tripped other studies in this area, which we discussed in our look at explainable AI for brain tumor MRI. Here the authors at least release code and source data, which makes an outside check possible.

What a practitioner can reuse

Three ideas travel well beyond this paper. The first is the labeling loop. If you have any task where a fluorescent or otherwise expensive channel reveals the truth, you can label with that channel and still deploy a model on the cheap one. The second is reading kinetics. A curve over time carries information that a single frame cannot, and it can tolerate a noisy per frame detector. The third is the small debris classifier. When detection errors have a simple geometric signature, a two feature model can clean them up at almost no cost.

The same ideas come with checks you should add. Define the peak rule numerically and publish the threshold. Report the spread across folds. Hold out entire compound classes, not just wells, when you claim generalization. And show per compound statistics for every DAMP, including the ones that fail to reach significance.

Conclusion

The achievement is concrete. Kim, Jang, Kim and colleagues, with Jong-Ho Cha and Hyungyu Lee as corresponding authors, built a screen that reads cell death as a story told over time and not as a single number. A detector trained with fluorescence and deployed on plain DIC images can count swelling cells and dead cells, and the resulting curves separated four apoptosis inducers from four swelling type killers in a blind panel.

The conceptual shift is the more interesting part. Earlier ICD screens measured one DAMP at one endpoint. This one watches the sequence unfold and asks whether swelling comes before rupture. That is closer to how ICD is described biologically, and it turns morphology into a screening signal that costs no reagents once the model exists.

The approach should transfer. Any process with a visible lifecycle and a fluorescent reporter for the training phase could use the same recipe, whether that is other regulated cell death forms, drug induced stress responses or differentiation. The authors suggest training the system to separate several kinds of cell death and to estimate how mixed populations divide among them. That would be valuable, since a single drug can trigger different death modes depending on context.

The remaining limitations are real. Swelling cell detection is noisy, with precision of 0.434 and mAP50 of 0.349. The blind test overlaps with the training compound and is dominated by ferroptosis inducers. Extracellular ATP supports the DAMP story for only one of four hits, and dosing is inconsistent across assays. These are fixable problems, and the open code makes fixing them possible for people outside the original lab.

The most useful next experiments are plain ones. Train on a second cell line and test on a third. Include ICD inducers from other mechanisms in a truly held out panel. Publish the peak calling rule and the fold by fold spread. If those checks hold, this becomes a practical front door for ICD screening. Until then it is a promising detector of swelling and rupture, and the gap between those two descriptions is where the next paper will be written.

A compact PyTorch version of the pipeline

The code below re implements the ideas in the paper in a form you can run on a laptop. It includes a small anchor free detector with a channel adapter for the merged and DIC only phases, GIoU plus BCE losses, the two stage fluorescence pretraining and DIC fine tuning routine, an ensemble that pools fold models, the MAL loop, precision, recall, F1 and mAP50 metrics, the boosted stump debris filter from Equation 1, and the curve reading step. The paper does not publish a loss, a backbone beyond naming YOLOv9, or a peak rule, so those parts are my own simplifications and are marked as such in the docstring. The smoke test trains on synthetic frames and checks every component. For the authors’ real implementation, use their GitHub repository.

"""
ICD-HTS pipeline (Kim et al., Computers in Biology and Medicine 196, 2025, 110727)
A compact PyTorch re-implementation of the ideas in the paper.

What is faithful to the paper
  * Two detection classes, swelling cell (SC) and dead cell (DC)
  * Fluorescence pretraining on merged DIC + RFP + GFP, then DIC only fine tuning (FPDF)
  * Model assisted labeling (MAL) loop with expert review
  * Ensemble of cross validation models
  * Debris filter driven by box width and height, Eq. (1), boosted trees with a sigmoid
  * Precision, recall, F1 and mAP50, and per frame SC / DC count curves

What is my own simplification (the paper does not publish these details)
  * A small anchor free YOLO style network instead of YOLOv9-e
  * GIoU box loss + BCE objectness + BCE class loss
  * Boosted decision stumps written in torch instead of calling XGBoost
  * The SC peak rule at the bottom, since the paper gives no numeric threshold
"""
import math
import random
from typing import Callable, List, Tuple

import torch
import torch.nn as nn
import torch.nn.functional as F

SC, DC = 0, 1
NUM_CLASSES = 2


# ----------------------------------------------------------------------------
# 1. Network building blocks
# ----------------------------------------------------------------------------
class ConvBNAct(nn.Module):
    def __init__(self, c_in, c_out, k=3, s=1):
        super().__init__()
        self.conv = nn.Conv2d(c_in, c_out, k, s, k // 2, bias=False)
        self.bn = nn.BatchNorm2d(c_out)
        self.act = nn.SiLU(inplace=True)

    def forward(self, x):
        return self.act(self.bn(self.conv(x)))


class CSPBlock(nn.Module):
    """Cross stage partial block, a light stand in for the YOLOv9 RepNCSP blocks."""

    def __init__(self, c_in, c_out, n=1):
        super().__init__()
        mid = c_out // 2
        self.split_a = ConvBNAct(c_in, mid, 1)
        self.split_b = ConvBNAct(c_in, mid, 1)
        self.blocks = nn.Sequential(*[ConvBNAct(mid, mid, 3) for _ in range(n)])
        self.merge = ConvBNAct(2 * mid, c_out, 1)

    def forward(self, x):
        return self.merge(torch.cat([self.blocks(self.split_a(x)), self.split_b(x)], 1))


class ChannelAdapter(nn.Module):
    """
    Feeds the same stem in both training phases.
    mode "merged" expects a 3 channel tensor (DIC, RFP, GFP) for pretraining.
    mode "dic" expects one channel and repeats it so pretrained stem weights are reused.
    """

    def __init__(self):
        super().__init__()
        self.mode = "merged"

    def forward(self, x):
        if self.mode == "dic" and x.shape[1] == 1:
            x = x.repeat(1, 3, 1, 1)
        return x


class ICDDetector(nn.Module):
    """Anchor free detector with one head at stride 8. Output per cell is
    (l, t, r, b distances in pixels, objectness logit, SC logit, DC logit)."""

    def __init__(self, width=24):
        super().__init__()
        self.adapter = ChannelAdapter()
        self.stem = ConvBNAct(3, width, 3, 2)
        self.stage1 = nn.Sequential(ConvBNAct(width, width * 2, 3, 2), CSPBlock(width * 2, width * 2))
        self.stage2 = nn.Sequential(ConvBNAct(width * 2, width * 4, 3, 2), CSPBlock(width * 4, width * 4, 2))
        self.neck = ConvBNAct(width * 4, width * 4, 3)
        self.box = nn.Conv2d(width * 4, 4, 1)
        self.obj = nn.Conv2d(width * 4, 1, 1)
        self.cls = nn.Conv2d(width * 4, NUM_CLASSES, 1)
        self.stride = 8
        nn.init.constant_(self.obj.bias, -4.0)  # start with few positives
        nn.init.constant_(self.cls.bias, -4.0)

    def set_mode(self, mode: str):
        assert mode in ("merged", "dic")
        self.adapter.mode = mode

    def forward(self, x):
        h = self.neck(self.stage2(self.stage1(self.stem(self.adapter(x)))))
        box = F.softplus(self.box(h)) * self.stride * 2  # positive distances
        return torch.cat([box, self.obj(h), self.cls(h)], 1)  # B x (4+1+C) x H x W


# ----------------------------------------------------------------------------
# 2. Geometry helpers, loss, decoding
# ----------------------------------------------------------------------------
def grid_centers(h, w, stride, device):
    ys, xs = torch.meshgrid(torch.arange(h, device=device), torch.arange(w, device=device), indexing="ij")
    return (xs.float() + 0.5) * stride, (ys.float() + 0.5) * stride


def box_iou(a, b):
    """a: N x 4, b: M x 4 in xyxy. Returns N x M IoU."""
    area_a = (a[:, 2] - a[:, 0]).clamp(min=0) * (a[:, 3] - a[:, 1]).clamp(min=0)
    area_b = (b[:, 2] - b[:, 0]).clamp(min=0) * (b[:, 3] - b[:, 1]).clamp(min=0)
    lt = torch.max(a[:, None, :2], b[None, :, :2])
    rb = torch.min(a[:, None, 2:], b[None, :, 2:])
    inter = (rb - lt).clamp(min=0).prod(-1)
    return inter / (area_a[:, None] + area_b[None] - inter + 1e-9)


def giou_loss(p, t):
    """Paired GIoU loss for N x 4 xyxy boxes."""
    lt = torch.max(p[:, :2], t[:, :2])
    rb = torch.min(p[:, 2:], t[:, 2:])
    inter = (rb - lt).clamp(min=0).prod(-1)
    ap = (p[:, 2] - p[:, 0]).clamp(min=0) * (p[:, 3] - p[:, 1]).clamp(min=0)
    at = (t[:, 2] - t[:, 0]) * (t[:, 3] - t[:, 1])
    union = ap + at - inter + 1e-9
    iou = inter / union
    clt = torch.min(p[:, :2], t[:, :2])
    crb = torch.max(p[:, 2:], t[:, 2:])
    carea = (crb - clt).prod(-1) + 1e-9
    return (1 - (iou - (carea - union) / carea)).mean()


def detection_loss(pred, targets: List[torch.Tensor], stride=8, w_box=5.0, w_obj=1.0, w_cls=1.0):
    """
    pred   B x 7 x H x W
    targets list of B tensors, each N x 5 with rows (cls, x1, y1, x2, y2) in pixels.
    A grid cell is positive when its centre falls inside a target box and is closest to that box centre.
    """
    B, _, H, W = pred.shape
    cx, cy = grid_centers(H, W, stride, pred.device)
    obj_t = torch.zeros(B, H, W, device=pred.device)
    cls_t = torch.zeros(B, NUM_CLASSES, H, W, device=pred.device)
    box_p, box_t = [], []
    for b in range(B):
        tg = targets[b]
        if tg.numel() == 0:
            continue
        x1, y1, x2, y2 = tg[:, 1], tg[:, 2], tg[:, 3], tg[:, 4]
        inside = (cx[None] > x1[:, None, None]) & (cx[None] < x2[:, None, None]) & \
                 (cy[None] > y1[:, None, None]) & (cy[None] < y2[:, None, None])
        tcx, tcy = (x1 + x2) / 2, (y1 + y2) / 2
        dist = (cx[None] - tcx[:, None, None]) ** 2 + (cy[None] - tcy[:, None, None]) ** 2
        dist = dist.masked_fill(~inside, float("inf"))
        best_d, best_i = dist.min(0)  # per cell, nearest containing box
        pos = torch.isfinite(best_d)
        if pos.sum() == 0:
            continue
        obj_t[b][pos] = 1.0
        ys, xs = pos.nonzero(as_tuple=True)
        ti = best_i[pos]
        cls_t[b][tg[ti, 0].long(), ys, xs] = 1.0
        d = pred[b, :4, ys, xs].T
        pcx, pcy = cx[ys, xs], cy[ys, xs]
        box_p.append(torch.stack([pcx - d[:, 0], pcy - d[:, 1], pcx + d[:, 2], pcy + d[:, 3]], 1))
        box_t.append(tg[ti, 1:5])
    l_obj = F.binary_cross_entropy_with_logits(pred[:, 4], obj_t, reduction="mean")
    pos_mask = obj_t.unsqueeze(1).expand_as(cls_t) > 0
    l_cls = F.binary_cross_entropy_with_logits(pred[:, 5:][pos_mask], cls_t[pos_mask]) if pos_mask.any() else pred.sum() * 0
    l_box = giou_loss(torch.cat(box_p), torch.cat(box_t)) if box_p else pred.sum() * 0
    return w_box * l_box + w_obj * l_obj + w_cls * l_cls


def nms(boxes, scores, thr=0.5):
    order = scores.argsort(descending=True)
    keep = []
    while order.numel() > 0:
        i = order[0]
        keep.append(i.item())
        if order.numel() == 1:
            break
        iou = box_iou(boxes[i][None], boxes[order[1:]])[0]
        order = order[1:][iou < thr]
    return torch.tensor(keep, dtype=torch.long)


@torch.no_grad()
def decode(pred, stride=8, conf=0.25, iou_thr=0.5):
    """Returns per image tensors with rows (cls, score, x1, y1, x2, y2)."""
    B, _, H, W = pred.shape
    cx, cy = grid_centers(H, W, stride, pred.device)
    out = []
    for b in range(B):
        d = pred[b, :4]
        boxes = torch.stack([cx - d[0], cy - d[1], cx + d[2], cy + d[3]], -1).reshape(-1, 4)
        obj = pred[b, 4].sigmoid().reshape(-1)
        cls = pred[b, 5:].sigmoid().reshape(NUM_CLASSES, -1)
        score, label = (cls * obj[None]).max(0)
        m = score > conf
        boxes, score, label = boxes[m], score[m], label[m]
        rows = []
        for c in range(NUM_CLASSES):
            idx = (label == c).nonzero(as_tuple=True)[0]
            if idx.numel() == 0:
                continue
            k = nms(boxes[idx], score[idx], iou_thr)
            sel = idx[k]
            rows.append(torch.cat([torch.full((len(sel), 1), float(c)), score[sel, None], boxes[sel]], 1))
        out.append(torch.cat(rows) if rows else torch.zeros(0, 6))
    return out


# ----------------------------------------------------------------------------
# 3. Metrics from Section 2.4.2  (precision, recall, F1, AP50, mAP50)
# ----------------------------------------------------------------------------
def match_detections(dets, gts, cls, iou_thr=0.5):
    """Greedy matching by score. Returns tp flags (sorted by score), scores, number of gts."""
    d = dets[dets[:, 0] == cls]
    g = gts[gts[:, 0] == cls][:, 1:5] if gts.numel() else torch.zeros(0, 4)
    d = d[d[:, 1].argsort(descending=True)]
    used = torch.zeros(len(g), dtype=torch.bool)
    tp = torch.zeros(len(d))
    if len(d) and len(g):
        iou = box_iou(d[:, 2:6], g)
        for i in range(len(d)):
            j = iou[i].masked_fill(used, -1).argmax()
            if iou[i, j] >= iou_thr:
                tp[i] = 1
                used[j] = True
    return tp, d[:, 1], len(g)


def evaluate(all_dets, all_gts, conf=0.25):
    """Returns dict with precision, recall, f1, ap50 for each class and mAP50."""
    res = {}
    for cls, name in ((SC, "SC"), (DC, "DC")):
        tps, scs, n_gt = [], [], 0
        for dets, gts in zip(all_dets, all_gts):
            tp, sc, n = match_detections(dets, gts, cls)
            tps.append(tp); scs.append(sc); n_gt += n
        tp = torch.cat(tps) if tps else torch.zeros(0)
        sc = torch.cat(scs) if scs else torch.zeros(0)
        order = sc.argsort(descending=True)
        tp = tp[order]
        ctp, cfp = tp.cumsum(0), (1 - tp).cumsum(0)
        rec = ctp / max(n_gt, 1)
        prec = ctp / (ctp + cfp).clamp(min=1)
        # all point interpolated area under the precision recall curve
        mrec = torch.cat([torch.zeros(1), rec, torch.ones(1)])
        mpre = torch.cat([torch.zeros(1), prec, torch.zeros(1)])
        for i in range(len(mpre) - 2, -1, -1):
            mpre[i] = max(mpre[i], mpre[i + 1])
        idx = (mrec[1:] != mrec[:-1]).nonzero(as_tuple=True)[0]
        ap = ((mrec[idx + 1] - mrec[idx]) * mpre[idx + 1]).sum().item()
        keep = sc[order] >= conf
        tp_c = tp[keep].sum().item()
        n_pred = int(keep.sum().item())
        p = tp_c / max(n_pred, 1)
        r = tp_c / max(n_gt, 1)
        res[name] = dict(precision=p, recall=r, f1=2 * p * r / max(p + r, 1e-9), ap50=ap)
    res["mAP50"] = (res["SC"]["ap50"] + res["DC"]["ap50"]) / 2
    return res


# ----------------------------------------------------------------------------
# 4. Training loop, FPDF (fluorescence pretraining, DIC fine tuning)
# ----------------------------------------------------------------------------
def train(model, images, targets, epochs=10, batch_size=4, lr=0.01, momentum=0.937, wd=5e-4, device="cpu"):
    """SGD settings follow Section 2.8 of the paper (lr 0.01, momentum 0.937, weight decay 0.0005)."""
    model.to(device).train()
    opt = torch.optim.SGD(model.parameters(), lr=lr, momentum=momentum, weight_decay=wd, nesterov=True)
    n = len(images)
    history = []
    for ep in range(epochs):
        order = list(range(n))
        random.shuffle(order)
        total = 0.0
        for s in range(0, n, batch_size):
            idx = order[s:s + batch_size]
            x = torch.stack([images[i] for i in idx]).to(device)
            t = [targets[i].to(device) for i in idx]
            loss = detection_loss(model(x), t)
            opt.zero_grad()
            loss.backward()
            torch.nn.utils.clip_grad_norm_(model.parameters(), 10.0)
            opt.step()
            total += loss.item() * len(idx)
        history.append(total / n)
    return history


def fpdf_train(model, merged_imgs, dic_imgs, targets, pre_epochs=10, ft_epochs=10):
    """Stage 1 on merged (DIC, RFP, GFP). Stage 2 on DIC only, keeping the learned weights."""
    model.set_mode("merged")
    h1 = train(model, merged_imgs, targets, epochs=pre_epochs)
    model.set_mode("dic")
    h2 = train(model, dic_imgs, targets, epochs=ft_epochs, lr=0.005)
    return h1, h2


@torch.no_grad()
def predict(model, images, conf=0.25, device="cpu"):
    model.to(device).eval()
    out = []
    for s in range(0, len(images), 8):
        out += decode(model(torch.stack(images[s:s + 8]).to(device)), conf=conf)
    return out


def ensemble_predict(models, images, conf=0.25, iou_thr=0.5):
    """Pools the boxes of the fold models and runs class wise NMS over the union."""
    per_model = [predict(m, images, conf) for m in models]
    fused = []
    for i in range(len(images)):
        d = torch.cat([pm[i] for pm in per_model])
        rows = []
        for c in (SC, DC):
            dc = d[d[:, 0] == c]
            if len(dc):
                rows.append(dc[nms(dc[:, 2:6], dc[:, 1], iou_thr)])
        fused.append(torch.cat(rows) if rows else torch.zeros(0, 6))
    return fused


# ----------------------------------------------------------------------------
# 5. Model assisted labeling (Fig. 3A)
# ----------------------------------------------------------------------------
def mal_round(model, unlabeled, annotator: Callable[[int, torch.Tensor], torch.Tensor], conf=0.3):
    """
    One MAL iteration. The model proposes pseudo labels, the annotator reviews each image and
    returns corrected labels (N x 5, cls x1 y1 x2 y2). Returns accepted labels and the
    number of proposals the expert had to remove or add, a proxy for review effort.
    """
    proposals = predict(model, unlabeled, conf)
    accepted, edits = [], 0
    for i, p in enumerate(proposals):
        pseudo = p[:, [0, 2, 3, 4, 5]]
        fixed = annotator(i, pseudo)
        edits += abs(len(fixed) - len(pseudo))
        accepted.append(fixed)
    return accepted, edits


# ----------------------------------------------------------------------------
# 6. Debris filter, Eq. (1): sigmoid( sum_k eta * f_k(x) ) >= theta  ->  SC else debris
# ----------------------------------------------------------------------------
class BoostedStumps:
    """Gradient boosted decision stumps on x = (width, height), logistic loss, Newton leaf values."""

    def __init__(self, n_trees=30, eta=0.3, reg_lambda=1.0, theta=0.5):
        self.n_trees, self.eta, self.lam, self.theta = n_trees, eta, reg_lambda, theta
        self.trees: List[Tuple[int, float, float, float]] = []  # feature, threshold, left, right
        self.base = 0.0

    def fit(self, X: torch.Tensor, y: torch.Tensor):
        p0 = y.float().mean().clamp(1e-3, 1 - 1e-3)
        self.base = torch.log(p0 / (1 - p0)).item()
        raw = torch.full_like(y, self.base, dtype=torch.float)
        for _ in range(self.n_trees):
            p = raw.sigmoid()
            g, h = p - y.float(), p * (1 - p)
            best = None
            for f in range(X.shape[1]):
                for thr in torch.unique(X[:, f]):
                    left = X[:, f] <= thr
                    if left.sum() == 0 or (~left).sum() == 0:
                        continue
                    gl, hl, gr, hr = g[left].sum(), h[left].sum(), g[~left].sum(), h[~left].sum()
                    gain = gl ** 2 / (hl + self.lam) + gr ** 2 / (hr + self.lam) - (gl + gr) ** 2 / (hl + hr + self.lam)
                    if best is None or gain > best[0]:
                        best = (gain.item(), f, thr.item(), (-gl / (hl + self.lam)).item(), (-gr / (hr + self.lam)).item())
            _, f, thr, vl, vr = best
            self.trees.append((f, thr, vl, vr))
            raw = raw + self.eta * torch.where(X[:, f] <= thr, torch.tensor(vl), torch.tensor(vr))
        return self

    def log_odds(self, X):
        raw = torch.full((len(X),), self.base)
        for f, thr, vl, vr in self.trees:
            raw = raw + self.eta * torch.where(X[:, f] <= thr, torch.tensor(vl), torch.tensor(vr))
        return raw

    def predict_proba(self, X):
        return self.log_odds(X).sigmoid()

    def is_sc(self, X):
        return self.predict_proba(X) >= self.theta


def filter_debris(dets: torch.Tensor, clf: BoostedStumps, image_size: float):
    """Drops SC detections the classifier calls debris. Width and height are normalised by image size."""
    if len(dets) == 0:
        return dets
    keep = torch.ones(len(dets), dtype=torch.bool)
    sc = dets[:, 0] == SC
    if sc.any():
        wh = torch.stack([dets[sc, 4] - dets[sc, 2], dets[sc, 5] - dets[sc, 3]], 1) / image_size
        keep[sc] = clf.is_sc(wh)
    return dets[keep]


# ----------------------------------------------------------------------------
# 7. From detections to the time series the paper reads (Fig. 5)
# ----------------------------------------------------------------------------
def count_curves(frame_dets: List[torch.Tensor]):
    sc = torch.tensor([(d[:, 0] == SC).sum().item() for d in frame_dets], dtype=torch.float)
    dc = torch.tensor([(d[:, 0] == DC).sum().item() for d in frame_dets], dtype=torch.float)
    return sc, dc


def has_sc_peak(sc_curve: torch.Tensor, dc_curve: torch.Tensor, min_height=3.0, ratio=0.15):
    """
    My own rule, because the paper does not state one. A compound is flagged when DC rises and the SC
    curve shows a transient peak, meaning its maximum is clearly above both its start and its end and
    is a meaningful fraction of the final DC count.
    """
    if dc_curve[-1] - dc_curve[0] < min_height:
        return False  # no death, so no call
    peak = sc_curve.max()
    shoulder = torch.max(sc_curve[:3].mean(), sc_curve[-3:].mean())
    return bool(peak - shoulder >= max(min_height, ratio * dc_curve[-1]))


# ----------------------------------------------------------------------------
# 8. Smoke test on synthetic data
# ----------------------------------------------------------------------------
def synth_frame(size=96, n_sc=2, n_dc=2, seed=0):
    """Fake DIC, RFP, GFP planes. SC are big discs with diffuse GFP, DC are small discs with RFP."""
    g = torch.Generator().manual_seed(seed)
    dic = 0.5 + 0.05 * torch.randn(size, size, generator=g)
    rfp = torch.zeros(size, size)
    gfp = torch.zeros(size, size)
    yy, xx = torch.meshgrid(torch.arange(size), torch.arange(size), indexing="ij")
    labels = []
    for cls, n, r in ((SC, n_sc, 11), (DC, n_dc, 7)):
        for _ in range(n):
            cx = int(torch.randint(r + 2, size - r - 2, (1,), generator=g))
            cy = int(torch.randint(r + 2, size - r - 2, (1,), generator=g))
            disc = ((xx - cx) ** 2 + (yy - cy) ** 2) <= r ** 2
            dic[disc] += 0.25 if cls == SC else -0.25
            (gfp if cls == SC else rfp)[disc] = 1.0
            labels.append([cls, cx - r, cy - r, cx + r, cy + r])
    return torch.stack([dic, rfp, gfp]), torch.tensor(labels, dtype=torch.float)


def smoke_test():
    torch.manual_seed(0)
    random.seed(0)
    size, n = 96, 24
    data = [synth_frame(size, 2, 2, seed=i) for i in range(n)]
    merged = [d[0] for d in data]
    dic = [d[0][:1] for d in data]
    labels = [d[1] for d in data]

    model = ICDDetector(width=16)
    out = model(torch.stack(merged[:2]))
    assert out.shape == (2, 4 + 1 + NUM_CLASSES, size // 8, size // 8), out.shape

    h1, h2 = fpdf_train(model, merged[:16], dic[:16], labels[:16], pre_epochs=25, ft_epochs=15)
    assert h2[-1] < h1[0], "loss should fall over the two stages"

    test_d, test_l = dic[16:], labels[16:]
    dets = predict(model, test_d, conf=0.2)
    metrics = evaluate(dets, test_l, conf=0.2)
    print("loss first/last pretrain", round(h1[0], 3), round(h1[-1], 3), "| finetune last", round(h2[-1], 3))
    print("SC", {k: round(v, 3) for k, v in metrics["SC"].items()})
    print("DC", {k: round(v, 3) for k, v in metrics["DC"].items()})
    print("mAP50", round(metrics["mAP50"], 3))

    ens = ensemble_predict([model, model], test_d, conf=0.2)
    assert all(len(e) <= len(d) for e, d in zip(ens, [torch.cat([x, x]) for x in dets]))

    # MAL round with a perfect annotator, only checks that the plumbing works
    gt = {i: labels[16 + i] for i in range(len(test_d))}
    _, edits = mal_round(model, test_d, lambda i, pseudo: gt[i], conf=0.2)
    print("MAL review edits on 8 images", edits)

    # Debris filter. SC are wide, debris are small.
    torch.manual_seed(1)
    sc_wh = 0.23 + 0.03 * torch.randn(200, 2)
    deb_wh = 0.06 + 0.02 * torch.randn(200, 2)
    X = torch.cat([sc_wh, deb_wh]).clamp(min=0.01)
    y = torch.cat([torch.ones(200), torch.zeros(200)])
    clf = BoostedStumps(n_trees=15).fit(X, y)
    acc = (clf.is_sc(X).float() == y).float().mean().item()
    print("debris filter accuracy", round(acc, 3))
    assert acc > 0.95
    fake = torch.tensor([[SC, 0.9, 10, 10, 32, 32], [SC, 0.8, 50, 50, 56, 56], [DC, 0.7, 60, 10, 74, 24]])
    kept = filter_debris(fake, clf, image_size=96.0)
    assert len(kept) == 2 and (kept[:, 0] == DC).sum() == 1  # the tiny SC is dropped, DC untouched

    # Time series rule on two fake compounds
    t = torch.arange(100, dtype=torch.float)
    dc_rise = 150 / (1 + torch.exp(-(t - 40) / 4))
    sc_peak = 25 * torch.exp(-((t - 38) ** 2) / 30)
    assert has_sc_peak(sc_peak, dc_rise) is True
    assert has_sc_peak(torch.zeros(100) + 1, dc_rise) is False       # death without swelling
    assert has_sc_peak(sc_peak, torch.zeros(100) + 2) is False       # nothing happens
    print("smoke test passed")


if __name__ == "__main__":
    smoke_test()

Frequently asked questions

Related reading