Neutrophil Reorientation via Local Feedback and PIP3-Driven Delay
Persistent Identifier
Use this permanent link to cite or share this Morpheus model:
How does immune cell polarity respond to changing stimuli?
Introduction
Neutrophils are fast-moving immune cells that polarize and reorient in response to external signals. Using optogenetic control of PI3K, recent experiments revealed that when a localized light stimulus becomes global, cells reverse their polarity — a behaviour not predicted by classical models.
This model suite tests a series of minimal feedback motifs that reproduce these dynamics using Cellular Potts Models (CPM) coupled to reaction-diffusion systems. The models include increasing layers of complexity: a base wave-pinning (WP) model, a WP model with a Rac-dependent inhibitor (WPI), and a full WPI-PIP3 model that captures the full experimental behaviour, including reversal.
Model Overview
We model the cell perimeter as a 1D membrane where active Rac ($u$), inactive Rac ($v$), a Rac-dependent inhibitor ($h$), and a PIP3-like activator ($p$) evolve under reaction-diffusion dynamics.
To simulate stimulus saturation, the light field $S(x,t)$ is normalized by its integral:
$$ \text{Effective stimulus} = \frac{S(x,t)}{\int S(x,t), dx} $$
This ensures that spatially broader stimuli do not lead to proportionally larger Rac activation, mimicking receptor saturation effects.
1. Wave-Pinning Model (WP)
Based on Mori et al. (2008), the classic WP model describes bistable Rac dynamics with conservation of total Rac:
$$ \frac{\partial u}{\partial t} = v \left( k_0 + \alpha \frac{S(x,t)}{\int S(x,t), dx} + \gamma \frac{u^n}{K^n + u^n} \right) - \delta u + D_u \nabla^2 u $$
$$ \frac{\partial v}{\partial t} = -v \left( k_0 + \alpha \frac{S(x,t)}{\int S(x,t), dx} + \gamma \frac{u^n}{K^n + u^n} \right) + \delta u + D_v \nabla^2 v $$
2. Wave-Pinning with Inhibitor (WPI)
We extend the WP model by introducing a Rac-dependent inhibitor $h(x,t)$ that diffuses slowly and suppresses Rac locally:
$$ \frac{\partial h}{\partial t} = k_h u - \delta_h h + D_h \nabla^2 h $$
The Rac equations are modified to include this feedback:
$$ \frac{\partial u}{\partial t} = v \left( k_0 + \alpha \frac{S(x,t)}{\int S(x,t), dx} + \gamma \frac{u^n}{K^n + u^n} \right) - (\delta + k_1 h) u + D_u \nabla^2 u $$
$$ \frac{\partial v}{\partial t} = -v \left( k_0 + \alpha \frac{S(x,t)}{\int S(x,t), dx} + \gamma \frac{u^n}{K^n + u^n} \right) + (\delta + k_1 h) u + D_v \nabla^2 v $$
3. Wave-Pinning with Inhibitor and PIP3 (WPI-PIP3)
To capture stimulus memory, we introduce a dynamic PIP3 signal $p(x,t)$ that activates Rac indirectly:
$$ \frac{\partial p}{\partial t} = (T_p - p)(k_p + \beta \frac{S(x,t)}{\int S(x,t), dx}) - \delta_p p + D_p \nabla^2 p $$
Rac activation now depends on $p(x,t)$:
$$ \frac{\partial u}{\partial t} = v \left( k_0 + \alpha p + \gamma \frac{u^n}{K^n + u^n} \right) - (\delta + k_1 h) u + D_u \nabla^2 u $$
$$ \frac{\partial v}{\partial t} = -v \left( k_0 + \alpha p + \gamma \frac{u^n}{K^n + u^n} \right) + (\delta + k_1 h) u + D_v \nabla^2 v $$
Model Files Included
WPI-PIP3_main.xml: Full model with Rac, inhibitor, and PIP3WPI.xml: Intermediate model with Rac and inhibitorWP.xml: Baseline wave-pinning model
Reference
This model is the original used in the publication, up to technical updates:
J. Algorta, J.P. Town, O.D. Weiner, L. Edelstein-Keshet: Investigating local negative feedback of Rac activity by mathematical models and cell motility simulations. bioRxiv, 2025.
The model originates from a preprint and thus has not been formally peer-reviewed.
Model
WPI-PIP3_main.xml
XML Preview
<?xml version='1.0' encoding='UTF-8'?>
<MorpheusModel version="4">
<Description>
<Title>Wave-Pinning with Inhibitor and PIP3-Mediated Input (WPI-PIP3)</Title>
<Details>
MorpheusModelID: https://identifiers.org/morpheus/M2073
Software: Morpheus (open source). Download from: https://morpheus.gitlab.io
Full title: Wave-Pinning Model with Inhibitor and PIP3-Mediated Input
Authors: J. Algorta
Submitters: J. Algorta
Date: 2025-05-15
Reference:
[1] J. Algorta, J.P. Town, O.D. Weiner, L. Edelstein-Keshet. "Investigating local negative feedback of Rac activity by mathematical models and cell motility simulations." bioRxiv, 2025. [Link](https://www.biorxiv.org/content/10.1101/2025.05.05.651928)
Comment:
This final model introduces a dynamic upstream field representing PIP3 (\( p(x,t) \)), which indirectly drives Rac activation. PIP3 introduces a time delay and memory to the system, allowing cells to reverse direction under local-to-global stimulus switching. This model captures all observed experimental behaviours.
</Details>
</Description>
<Global>
<Constant symbol="ts1" name="Time to *end* first light stimulus from "time = 0"" value="0"/>
<Constant symbol="ts2" name="Time to *starts* second light stimulus from "time = ts2"" value="60"/>
<Constant symbol="L" name="Angular size of the light stimulus" value="pi/4"/>
<Variable symbol="A" value="0.0"/>
<Variable symbol="B" value="0.0"/>
<Variable symbol="U" value="-1"/>
<Field symbol="prot" value="0.0"/>
<Variable symbol="ts3" value="3000"/>
<Constant symbol="middle_temp" value="rand_uni(pi/4,2*pi-pi/4)"/>
<Constant symbol="lengthscale" value="2*pi/(2*pi*sqrt(2000/pi))"/>
<Variable symbol="vol_target" value="2000"/>
<!-- <Disabled>
<Variable symbol="vol_target" value="rand_norm(2000.142, 683.505) "/>
</Disabled>
-->
<Constant symbol="middle" value="3*pi/2"/>
<!-- <Disabled>
<Constant symbol="middle" value="middle_temp"/>
</Disabled>
-->
</Global>
<Space>
<Lattice class="square">
<Size symbol="lattice" value="500, 500, 0"/>
<BoundaryConditions>
<Condition type="periodic" boundary="x"/>
<Condition type="periodic" boundary="y"/>
</BoundaryConditions>
<NodeLength value="1.0"/>
<Neighborhood>
<Distance>2.0</Distance>
</Neighborhood>
</Lattice>
<SpaceSymbol symbol="l"/>
<MembraneLattice>
<Resolution symbol="memsize" value="100"/>
<SpaceSymbol symbol="m"/>
</MembraneLattice>
</Space>
<Time>
<StartTime value="0"/>
<StopTime value="400"/>
<SaveInterval value="0"/>
<TimeSymbol symbol="time"/>
<!-- <Disabled>
<RandomSeed value="1743554059"/>
</Disabled>
-->
</Time>
<CellTypes>
<CellType name="cells" class="biological">
<VolumeConstraint target="vol_target" strength="1"/>
<SurfaceConstraint target="1.46" strength="1" mode="aspherity"/>
<MembraneProperty symbol="A" name="Active form of protein" value="2.05*(m.phi>middle-pi/4)*(m.phi<middle+pi/4)">
<Diffusion rate="0.080/lengthscale^2"/>
</MembraneProperty>
<MembraneProperty symbol="B" name="Inactive form of protein" value="7.31-2/4">
<Diffusion rate="0.51/lengthscale^2"/>
</MembraneProperty>
<MembraneProperty symbol="H" name="Local Rac inhibtor" value="2">
<Diffusion rate=" 0.0374/lengthscale^2"/>
</MembraneProperty>
<MembraneProperty symbol="U" name="Light stimulus around the cell's membrane" value="0">
<Diffusion rate="0.0"/>
</MembraneProperty>
<MembraneProperty symbol="P3" name="PIP3" value="0.93">
<Diffusion rate="0.62/lengthscale^2"/>
</MembraneProperty>
<System time-scaling="2" name="WavePinning with Inhibitor" solver="Euler [fixed, O(1)]">
<Constant symbol="k_0" value="rand_norm(0.028, 0.000765)"/>
<Constant symbol="k_1" value="rand_norm(0.372, 0.006632)"/>
<Constant symbol="s_1" value="rand_norm(0.059, 0.001530)"/>
<Constant symbol="gamma" value="rand_norm(1.38, 0.0214)"/>
<Constant symbol="delta" value="rand_norm(0.721, 0.006887)"/>
<Constant symbol="K" value="rand_norm(2.93, 0.0408)"/>
<Constant symbol="n" name="Hill coefficient" value="2"/>
<Intermediate symbol="F" value="1*(B*(k_0+ s_1*(P3) + (gamma*A^n) / (K^n + A^n) ) - (delta+k_1*H)*A)"/>
<DiffEqn symbol-ref="A">
<Expression>F</Expression>
</DiffEqn>
<DiffEqn symbol-ref="B">
<Expression>-F</Expression>
</DiffEqn>
</System>
<System time-scaling="2" name="Inhibitor" solver="Euler [fixed, O(1)]">
<Constant symbol="alpha" value="rand_norm(0.0068, 0.000102040)"/>
<Constant symbol="deltah" value="rand_norm(0.004, 0.00015306)"/>
<DiffEqn symbol-ref="H">
<Expression>1*(alpha*A-deltah*H)</Expression>
</DiffEqn>
</System>
<System time-scaling="0.004" name="PIP3" solver="Euler [fixed, O(1)]">
<Intermediate symbol="u_damp" value="max(0.001,u_inte)"/>
<Intermediate symbol="G" value="((3.08-P3) * (0.41 + 0.26 * 1*3.7*U/(u_damp)+0.26) - (1.59) * P3)"/>
<DiffEqn symbol-ref="P3">
<Expression>G</Expression>
</DiffEqn>
</System>
<System name="Input system" solver="Dormand-Prince [adaptive, O(5)]">
<Annotation>Below are the different ways to chaneg how the light stimulus (U) is applied on the cell's membrane. ONLY ONE SHOULD BE TURNED ON AT A TIME
</Annotation>
<Rule symbol-ref="U" name="Local to Global">
<Expression>((abs(m.phi-3*pi/2)<L)*(time<ts1))+((abs(m.phi-3*pi/2+0.4)<5000)*(time>ts3))+((abs(m.phi-(p.phi+pi/2))<L)*(time>ts2)*(time<=ts3))*((p.phi+pi/2)<(2*pi-L))+(((abs(m.phi-(p.phi+pi/2))<L)*(time>ts2)*(time<=ts3))+((abs(m.phi-(p.phi-2*pi+pi/2))<L)*(time>ts2)*(time<=ts3)))*((p.phi+pi/2)>(2*pi-L))</Expression>
</Rule>
<!-- <Disabled>
<Rule symbol-ref="U" name="Side-Side">
<Expression>(((abs(m.phi - (angle + pi/2)) < L) * ((angle + pi/2) <= 2*pi - L)) +
((abs(m.phi - (angle + pi/2 - 2*pi)) < L) * ((angle + pi/2) > 2*pi - L) * ((angle + pi/2) < 2*pi + L)))*(time > ts2) * (time < ts3)
</Expression>
</Rule>
</Disabled>
-->
</System>
<PropertyVector symbol="p" name="polarity" value="0.0, 0.0, 0.0"/>
<MotilityReporter time-step="10" name="Velocity extraction">
<Velocity symbol-ref="p"/>
</MotilityReporter>
<StarConvex membrane="A" strength="30" name="Protrusive domain identification"/>
<Protrusion strength="10" field="prot" name="Protrusion persistence" maximum="2500"/>
<Event name="Trigger condition">
<Condition>(abs(p.phi-pi/2)<0.3)*(time>ts2+10)</Condition>
<Rule symbol-ref="ts3">
<Expression>time+0</Expression>
</Rule>
</Event>
<Mapper name="Light direction">
<Input value="U"/>
<Polarity symbol-ref="p1"/>
</Mapper>
<PropertyVector symbol="p1" name="polarity1" value="0.0, 0.0, 0.0"/>
<ConnectivityConstraint name=""/>
<Property symbol="u_inte" name="Intermediate variable" value="0.0"/>
<Mapper name="Integral of Light">
<Input value="2*pi*U/memsize"/>
<Output symbol-ref="u_inte" mapping="sum"/>
</Mapper>
</CellType>
</CellTypes>
<CPM>
<Interaction default="0.0"/>
<MonteCarloSampler stepper="edgelist">
<MCSDuration value="0.07"/>
<Neighborhood>
<Order>2</Order>
</Neighborhood>
<MetropolisKinetics yield="0.05" temperature="0.3"/>
</MonteCarloSampler>
<ShapeSurface scaling="norm">
<Neighborhood>
<Order>6</Order>
</Neighborhood>
</ShapeSurface>
</CPM>
<CellPopulations>
<Population type="cells" size="0">
<InitCellObjects mode="distance">
<Arrangement repetitions="1, 1, 0" displacements="1, 1, 0">
<Sphere center="100 250 0" radius="24"/>
</Arrangement>
</InitCellObjects>
</Population>
</CellPopulations>
<Analysis>
<Gnuplotter time-step="5" decorate="false">
<Terminal name="png"/>
<Plot title="Active Protein & Light Stimulus">
<Cells flooding="true" value="A"/>
<CellArrows orientation="p*10"/>
<Field symbol-ref="U" max="1" min="0">
<ColorMap>
<Color color="white" value="0"/>
<Color color="blue" value="1"/>
</ColorMap>
</Field>
</Plot>
</Gnuplotter>
<CellTracker time-step="1" format="ISBI 2012 (XML)"/>
<Logger time-step="1">
<Input>
<Symbol symbol-ref="m.phi"/>
<Symbol symbol-ref="A"/>
<Symbol symbol-ref="B"/>
<Symbol symbol-ref="U"/>
<Symbol symbol-ref="p.phi"/>
<Symbol symbol-ref="H"/>
<Symbol symbol-ref="P3"/>
</Input>
<Output>
<TextOutput/>
</Output>
<Plots>
<Plot time-step="5">
<Style style="points"/>
<Terminal terminal="png"/>
<X-axis maximum="6.283">
<Symbol symbol-ref="m.phi"/>
</X-axis>
<Y-axis minimum="0.0" maximum="7">
<Symbol symbol-ref="B"/>
<Symbol symbol-ref="A"/>
<Symbol symbol-ref="U"/>
<Symbol symbol-ref="H"/>
<Symbol symbol-ref="P3"/>
</Y-axis>
<Range>
<Time mode="current"/>
</Range>
</Plot>
</Plots>
</Logger>
<ModelGraph include-tags="#untagged" format="dot" reduced="false"/>
<!-- <Disabled>
<Logger time-step="1.0">
<Input>
<Symbol symbol-ref="ts3"/>
<Symbol symbol-ref="p.phi"/>
</Input>
<Output>
<TextOutput/>
</Output>
</Logger>
</Disabled>
-->
</Analysis>
</MorpheusModel>
Downloads
Files associated with this model: