Skip to content

virus antibody model #253

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 7 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
61 changes: 61 additions & 0 deletions examples/virus_antibody/README.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
# Virus-Antibody Model

This model is a simulation of immune reaction declined as a confrontation between antibody agents and virus agents. The global idea is to model how the immune system can struggle against new virus but is able to adapt over time and beat a same virus if it comes back. The results are quite interesting as the simulation can go both ways (virus win or antibodies win) with a little tweak in the base parameters.


**It showcases :**
- **Usage of memory in agents** : divided into a short term memory using a deque to easily add and remove memories in case of a new virus encounter, and a long term memory (here a simple list)
- **Agent knowledge sharing** : the antibodies are able to share short term memory)
- **Usage of weak referencing** to avoid coding errors (antibodies can store viruses in a `self.target` attribute)
- Emergence of completely **different outcomes** with only small changes in parameters


For example, with a given set of fixed parameters :
| Virus mutation rate = 0.15 (antibodies win) | Virus mutation rate = 0.2 (viruses win) |
|--------------------------------------------------|--------------------------------------------------|
| ![](images/antibodies_win.png) | ![](images/viruses_win.png) |




## How It Works

1. **Initialization**: The model initializes a population of viruses and antibodies in a continuous 2D space.
2. **Agent Behavior**:
- Antibodies move randomly until they detect a virus within their sight range (becomes purple), than pursue the virus.
- Antibodies pass on all the virus DNA in their short term memory to the nearest antibodies (cf. example)
- Viruses move randomly and can duplicate or mutate.
3. **Engagement (antibody vs virus)**: When an antibody encounters a virus:
- If the antibody has the virus's DNA in its memory, it destroys the virus.
- Otherwise, the virus may defeat the antibody, causing it to lose health or become inactive temporarily.
4. **Duplication**: Antibodies and viruses can duplicate according to their duplication rate.


> Example for memory transmission : Let's look at two antibodies A1 and A2
> `A1.st_memory() = [ABC]` and `A1.lt_memory() = [ABC]`
> `A2.st_memory() = [DEF]` and `A2.lt() = [DEF]`
>
> After A1 encounters A2,
> `A1.st_memory() = [DEF]` and `A1.lt() = [ABC, DEF]`
> `A2.st_memory() = [ABC]` and `A1.lt() = [DEF, ABC]`
>
> A1 and A2 'switched' short term memory but both have the two viruses DNA in their long term memory

For further details, here is the full architecture of this model :

<div align="center">
<img src="images/virus_antibody_architecture.png" width="550"/>
</div>

## Usage

After cloning the repo and installing mesa on pip, run the application with :
```bash
solara run app.py
```

## A couple more of interesting cases

| An interesting tendency inversion | high duplication + high mutation = both grow (more viruses) | high duplication + low mutation = both grow (more antibodies) |
|---|---|---|
| <img src="images/pattern.png" width="550"/> | <img src="images/grow_virus_wins.png" width="450"/> | <img src="images/grow_antibody_wins.png" width="450"/> |
259 changes: 259 additions & 0 deletions examples/virus_antibody/agents.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,259 @@
"""
Mesa implementation of Virus/Antibody model: Agents module.
"""

import copy
import os
import sys
import weakref
from collections import deque

import numpy as np

sys.path.insert(0, os.path.abspath("../../mesa"))
from mesa.experimental.continuous_space import ContinuousSpaceAgent


class AntibodyAgent(ContinuousSpaceAgent):
"""An Antibody agent. They move randomly until they see a virus, go fight it.
If they lose, stay KO for a bit, lose health and back to random moving.
"""

def __init__(
self,
model,
space,
sight_range,
duplication_rate,
ko_timeout,
memory_capacity,
initial_position=(0, 0),
direction=(1, 1),
):
super().__init__(model=model, space=space)

# Movement & state
self.position = initial_position
self.speed = 1.5
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it be possible to make speed and health kwargs on the init? In my view it is good practice to make all attributes explicitly controllable.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I made the choice of limiting the parameter to these ones because otherwise it's a lot for the user in my opinion. If I add all of them, it would mean more than 10 parameters to be decided by the user. Do you think I should still add all of the other parameters ? Or are there one or two more that would be good to add and still leave the others hidden (like ko_timeout, sight_range, and a few non-crucial parameters) ?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fair enough. One option is to turn the constants that are the same for all agents of a given class into class level attributes. This still makes it easy to change them without having to add them all to the signature of the __init__.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good idea ! It does allow to delete quite a few lines of code in model.py also.

self.direction = np.array(direction, dtype=float)

# Characteristics
self.sight_range = sight_range
self.health = 2
self.duplication_rate = duplication_rate

# Memory
self.st_memory: deque = deque()
self.lt_memory: list = []
self.memory_capacity = memory_capacity

# Target & KO state
self.target = None # will hold a weakref.ref or None
self.ko_timeout = ko_timeout
self.ko_steps_left = 0

def step(self):
# Acquire a virus target if we don't already have one
if self.target is None:
closest = self.find_closest_virus()
if closest:
self.target = weakref.ref(closest)

# Communicate and maybe duplicate
self.communicate()
if self.random.random() < self.duplication_rate:
self.duplicate()

# Then move
self.move()

def find_closest_virus(self):
agents, _ = self.space.get_agents_in_radius(self.position, self.sight_range)
viruses = [a for a in agents if isinstance(a, VirusAgent)]
return viruses[0] if viruses else None

def communicate(self) -> bool:
agents, _ = self.space.get_agents_in_radius(self.position, self.sight_range)
peers = [
a
for a in agents
if isinstance(a, AntibodyAgent) and a.unique_id != self.unique_id
]
if not peers:
return False

for other in peers:
to_share = [
dna for dna in self.st_memory if dna and dna not in other.lt_memory
]
if to_share:
other.st_memory.extend(to_share)
other.lt_memory.extend(to_share)
while len(other.st_memory) > self.memory_capacity:
other.st_memory.popleft()
return True

def duplicate(self):
clone = AntibodyAgent(
self.model,
self.space,
sight_range=self.sight_range,
duplication_rate=self.duplication_rate,
ko_timeout=self.ko_timeout,
memory_capacity=self.memory_capacity,
initial_position=self.position,
direction=self.direction,
)
# Copy over memory
clone.st_memory = deque(item for item in self.st_memory if item)
clone.lt_memory = [item for item in self.lt_memory if item]
clone.target = None
clone.ko_steps_left = 0

self.model.antibodies_set.add(clone)

def move(self):
# If we've been removed from the space, bail out
if getattr(self, "space", None) is None:
return

# Dereference weakref if needed
target = (
self.target()
if isinstance(self.target, weakref.ReferenceType)
else self.target
)

new_pos = None

# KO state: target refers back to self
if target is self:
self.ko_steps_left -= 1
if self.ko_steps_left <= 0:
self.target = None

# Random walk if no target
elif target is None:
perturb = np.array(
[
self.random.uniform(-0.5, 0.5),
self.random.uniform(-0.5, 0.5),
]
)
self.direction = self.direction + perturb
norm = np.linalg.norm(self.direction)
if norm > 0:
self.direction /= norm
new_pos = self.position + self.direction * self.speed

# Chase a valid virus target
else:
if getattr(target, "space", None) is not None:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

again, this should not be needed if everything else is implemented correctly.

vec = np.array(target.position) - np.array(self.position)
dist = np.linalg.norm(vec)
if dist > self.speed:
self.direction = vec / dist
new_pos = self.position + self.direction * self.speed
else:
self.engage_virus(target)
else:
self.target = None

if new_pos is not None:
self.position = new_pos

def engage_virus(self, virus) -> str:
# If it's already gone
if virus not in self.model.agents:
self.target = None
return "no_target"

dna = copy.deepcopy(virus.dna)
if dna in self.st_memory or dna in self.lt_memory:
virus.remove()
self.target = None
return "win"
else:
# KO (or death)
self.health -= 1
if self.health <= 0:
self.remove()
return "dead"

self.st_memory.append(dna)
self.lt_memory.append(dna)
self.ko_steps_left = self.ko_timeout
# mark KO state by weak-ref back to self
self.target = weakref.ref(self)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I dont understand what is going on here.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When an antibody looses a confrontation with a virus (because the virus DNA is not in it's memory yet), the antibody is ko for a few steps. It can't move during this time and I symbolise this by setting the agent's target to itself.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, clear.

This might be implemented using event scheduling, which would result in a modest speedup of runtime.

return "ko"


class VirusAgent(ContinuousSpaceAgent):
"""A virus agent: random movement, mutation, duplication, passive to antibodies."""

def __init__(
self,
model,
space,
mutation_rate,
duplication_rate,
position=(0, 0),
dna=None,
):
super().__init__(model=model, space=space)

self.position = position
self.mutation_rate = mutation_rate
self.duplication_rate = duplication_rate
self.speed = 1
self.direction = np.array((1, 1), dtype=float)
self.dna = dna if dna is not None else self.generate_dna()

def step(self):
# If already removed from the space, don't do anything
if getattr(self, "space", None) is None:
return
if self.random.random() < self.duplication_rate:
self.duplicate()
self.move()

def duplicate(self):
clone = VirusAgent(
self.model,
self.space,
mutation_rate=self.mutation_rate,
duplication_rate=self.duplication_rate,
position=self.position,
dna=self.generate_dna(self.dna),
)
self.model.viruses_set.add(clone)

def generate_dna(self, dna=None):
if dna is None:
return [self.random.randint(0, 9) for _ in range(3)]
idx = self.random.randint(0, 2)
chance = self.random.random()
if chance < self.mutation_rate / 2:
dna[idx] = (dna[idx] + 1) % 10
elif chance < self.mutation_rate:
dna[idx] = (dna[idx] - 1) % 10
return dna

def move(self):
if getattr(self, "space", None) is None:
return

# Random walk
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Both agents have a random walk component. Might it be possible to implement this the same way for both and move it to a new super agent class from which both current agents derive?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It could be a good idea ! At first I thought it could also be used for the duplicate() method. However, there are a lot of parameters to take into account and the code is very short inside, so I don't think that it would really help to make it derive from a parent class.

So I the end, there is just a few lines for the random movement. Isn't a whole parent class a bit overkill ?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the choice is between repeated code or a parent class, I prefer the parent class. Of course, in python you could also move the relevant code into a helper function (not really appropriate for random walks).

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Alright, will do, thank you for the advice !
On another note, did you have time to check out why step was being called on a removed agent ?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Monday was a banking holiday so I had family obligations. I'll try to get to it today, but no promises unfortunately.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Of course, nothing urgent ! Only if you have time :)

perturb = np.array(
[
self.random.uniform(-0.5, 0.5),
self.random.uniform(-0.5, 0.5),
]
)
self.direction = self.direction + perturb
norm = np.linalg.norm(self.direction)
if norm > 0:
self.direction /= norm

# Step
self.position = self.position + self.direction * self.speed
Loading
Loading