diff --git a/core/PhysiCell_cell.cpp b/core/PhysiCell_cell.cpp
index 1f3e2646b..2c65c2190 100644
--- a/core/PhysiCell_cell.cpp
+++ b/core/PhysiCell_cell.cpp
@@ -238,12 +238,15 @@ Cell_State::Cell_State()
 	neighbors.resize(0); 
 	spring_attachments.resize(0); 
 
+	attacked_by.resize(0);
+
 	orientation.resize( 3 , 0.0 ); 
 	
 	simple_pressure = 0.0; 
 	
 	attached_cells.clear(); 
 	spring_attachments.clear(); 
+	attacked_by.clear();
 	
 	number_of_nuclei = 1; 
 	
@@ -457,6 +460,8 @@ Cell::~Cell()
 		{
 			// release any attached cells (as of 1.7.2 release)
 			this->remove_all_attached_cells(); 
+			this->remove_all_attackers();
+			this->remove_self_from_attacked();
 			// 1.11.0
 			this->remove_all_spring_attachments(); 
 
@@ -549,8 +554,12 @@ Cell* Cell::divide( )
 	// phenotype.flagged_for_division = false; 
 	// phenotype.flagged_for_removal = false; 
 	
-	// make sure ot remove adhesions 
+	// make sure to remove adhesions 
 	remove_all_attached_cells(); 
+	// the springs below include the one holding an attack together, so end the
+	// attack rather than leave it spring-less
+	remove_all_attackers();
+	remove_self_from_attacked();
 	remove_all_spring_attachments(); 
 
 	// version 1.10.3: 
@@ -651,6 +660,7 @@ Cell* Cell::divide( )
 	set_total_volume(phenotype.volume.total);
 	
 	// child->set_phenotype( phenotype ); 
+	// pAttackTarget was cleared above, so the daughter inherits a NULL one
 	child->phenotype = phenotype; 
 
     if (child->phenotype.intracellular){
@@ -1144,6 +1154,14 @@ void Cell::convert_to_cell_definition( Cell_Definition& cd )
 	Geometry cell_geometry = phenotype.geometry;
 	Molecular cell_molecular = phenotype.molecular;
 	Custom_Cell_Data cell_custom_data = custom_data;
+
+	// pAttackTarget is phenotype state, and the assignment below replaces it, so
+	// end this cell's own attack. attacked_by and the spring are in state, which
+	// transformation preserves, so attacks against this cell continue.
+	remove_self_from_attacked();
+	// should we also remove all attackers?  That would be a change in behavior, so for now we don't.
+	// remove_all_attackers();
+
 	// use the cell defaults; 
 	type = cd.type; 
 	type_name = cd.name; 
@@ -1210,6 +1228,8 @@ void delete_cell( int index )
 	
 	// release any attached cells (as of 1.7.2 release)
 	pDeleteMe->remove_all_attached_cells(); 
+	pDeleteMe->remove_all_attackers();
+	pDeleteMe->remove_self_from_attacked();
 	// 1.11.0 
 	pDeleteMe->remove_all_spring_attachments(); 
 
@@ -1244,6 +1264,8 @@ void delete_cell_original( int index ) // before June 11, 2020
 	
 	// release any attached cells (as of 1.7.2 release)
 	(*all_cells)[index]->remove_all_attached_cells(); 
+	(*all_cells)[index]->remove_all_attackers();
+	(*all_cells)[index]->remove_self_from_attacked();
 	// 1.11.0
 	(*all_cells)[index]->remove_all_spring_attachments(); 
 	
@@ -1347,133 +1369,136 @@ void Cell::ingest_cell( Cell* pCell_to_eat )
 	if( pCell_to_eat == this )
 	{ return; } 
 	
-	// don't ingest a cell that's already ingested 
-	if( pCell_to_eat->phenotype.volume.total < 1e-15 )
-	{ return; } 
-		
 	// make this thread safe 
 	#pragma omp critical
 	{
-		/*
-		if( pCell_to_eat->phenotype.death.dead == true )
-		{ std::cout << this->type_name << " (" << this << ")" << " eats dead " << pCell_to_eat->type_name << " (" << pCell_to_eat 
-			<< ") of size " << pCell_to_eat->phenotype.volume.total << std::endl; }
-		else
-		{ std::cout << this->type_name << " (" << this << ")" << " eats live " << pCell_to_eat->type_name << " (" << pCell_to_eat 
-			<< ") of size " << pCell_to_eat->phenotype.volume.total << std::endl; }
-		*/
-
-		// mark it as dead 
-		pCell_to_eat->phenotype.death.dead = true; 
-		// set secretion and uptake to zero 
-		pCell_to_eat->phenotype.secretion.set_all_secretion_to_zero( );  
-		pCell_to_eat->phenotype.secretion.set_all_uptake_to_zero( ); 
-		
-		// deactivate all custom function 
-		pCell_to_eat->functions.custom_cell_rule = NULL; 
-		pCell_to_eat->functions.update_phenotype = NULL; 
-		pCell_to_eat->functions.contact_function = NULL; 
-		pCell_to_eat->functions.cell_division_function = NULL; 
-		pCell_to_eat->functions.cell_division_direction_function = NULL; 
-		
-		// should set volume fuction to NULL too! 
-		pCell_to_eat->functions.volume_update_function = NULL; 
+		// don't ingest a cell that's already ingested 
+		if( pCell_to_eat->phenotype.volume.total >= 1e-15 )
+		{ 
+			/*
+			if( pCell_to_eat->phenotype.death.dead == true )
+			{ std::cout << this->type_name << " (" << this << ")" << " eats dead " << pCell_to_eat->type_name << " (" << pCell_to_eat 
+				<< ") of size " << pCell_to_eat->phenotype.volume.total << std::endl; }
+			else
+			{ std::cout << this->type_name << " (" << this << ")" << " eats live " << pCell_to_eat->type_name << " (" << pCell_to_eat 
+				<< ") of size " << pCell_to_eat->phenotype.volume.total << std::endl; }
+			*/
+
+			// mark it as dead 
+			pCell_to_eat->phenotype.death.dead = true; 
+			// set secretion and uptake to zero 
+			pCell_to_eat->phenotype.secretion.set_all_secretion_to_zero( );  
+			pCell_to_eat->phenotype.secretion.set_all_uptake_to_zero( ); 
+			
+			// deactivate all custom function 
+			pCell_to_eat->functions.custom_cell_rule = NULL; 
+			pCell_to_eat->functions.update_phenotype = NULL; 
+			pCell_to_eat->functions.contact_function = NULL; 
+			pCell_to_eat->functions.cell_division_function = NULL; 
+			pCell_to_eat->functions.cell_division_direction_function = NULL; 
+			
+			// should set volume fuction to NULL too! 
+			pCell_to_eat->functions.volume_update_function = NULL; 
 
-		// set cell as unmovable and non-secreting 
-		pCell_to_eat->is_movable = false; 
-		pCell_to_eat->is_active = false; 
+			// set cell as unmovable and non-secreting 
+			pCell_to_eat->is_movable = false; 
+			pCell_to_eat->is_active = false; 
 
-		// absorb all the volume(s)
+			// absorb all the volume(s)
 
-		// absorb fluid volume (all into the cytoplasm) 
-		phenotype.volume.cytoplasmic_fluid += pCell_to_eat->phenotype.volume.fluid; 
-		pCell_to_eat->phenotype.volume.cytoplasmic_fluid = 0.0; 
-		
-		// absorb nuclear and cyto solid volume (into the cytoplasm) 
-		phenotype.volume.cytoplasmic_solid += pCell_to_eat->phenotype.volume.cytoplasmic_solid; 
-		pCell_to_eat->phenotype.volume.cytoplasmic_solid = 0.0; 
-		
-		phenotype.volume.cytoplasmic_solid += pCell_to_eat->phenotype.volume.nuclear_solid; 
-		pCell_to_eat->phenotype.volume.nuclear_solid = 0.0; 
-		
-		// consistency calculations 
-		
-		phenotype.volume.fluid = phenotype.volume.nuclear_fluid + 
-			phenotype.volume.cytoplasmic_fluid; 
-		pCell_to_eat->phenotype.volume.fluid = 0.0; 
-		
-		phenotype.volume.solid = phenotype.volume.cytoplasmic_solid + 
-			phenotype.volume.nuclear_solid; 
-		pCell_to_eat->phenotype.volume.solid = 0.0; 
-		
-		// no change to nuclear volume (initially) 
-		pCell_to_eat->phenotype.volume.nuclear = 0.0; 
-		pCell_to_eat->phenotype.volume.nuclear_fluid = 0.0; 
-		
-		phenotype.volume.cytoplasmic = phenotype.volume.cytoplasmic_solid + 
-			phenotype.volume.cytoplasmic_fluid; 
-		pCell_to_eat->phenotype.volume.cytoplasmic = 0.0; 
-		
-		phenotype.volume.total = phenotype.volume.nuclear + 
-			phenotype.volume.cytoplasmic; 
-		pCell_to_eat->phenotype.volume.total = 0.0; 
+			// absorb fluid volume (all into the cytoplasm) 
+			phenotype.volume.cytoplasmic_fluid += pCell_to_eat->phenotype.volume.fluid; 
+			pCell_to_eat->phenotype.volume.cytoplasmic_fluid = 0.0; 
+			
+			// absorb nuclear and cyto solid volume (into the cytoplasm) 
+			phenotype.volume.cytoplasmic_solid += pCell_to_eat->phenotype.volume.cytoplasmic_solid; 
+			pCell_to_eat->phenotype.volume.cytoplasmic_solid = 0.0; 
+			
+			phenotype.volume.cytoplasmic_solid += pCell_to_eat->phenotype.volume.nuclear_solid; 
+			pCell_to_eat->phenotype.volume.nuclear_solid = 0.0; 
+			
+			// consistency calculations 
+			
+			phenotype.volume.fluid = phenotype.volume.nuclear_fluid + 
+				phenotype.volume.cytoplasmic_fluid; 
+			pCell_to_eat->phenotype.volume.fluid = 0.0; 
+			
+			phenotype.volume.solid = phenotype.volume.cytoplasmic_solid + 
+				phenotype.volume.nuclear_solid; 
+			pCell_to_eat->phenotype.volume.solid = 0.0; 
+			
+			// no change to nuclear volume (initially) 
+			pCell_to_eat->phenotype.volume.nuclear = 0.0; 
+			pCell_to_eat->phenotype.volume.nuclear_fluid = 0.0; 
+			
+			phenotype.volume.cytoplasmic = phenotype.volume.cytoplasmic_solid + 
+				phenotype.volume.cytoplasmic_fluid; 
+			pCell_to_eat->phenotype.volume.cytoplasmic = 0.0; 
+			
+			phenotype.volume.total = phenotype.volume.nuclear + 
+				phenotype.volume.cytoplasmic; 
+			pCell_to_eat->phenotype.volume.total = 0.0; 
 
-		phenotype.volume.fluid_fraction = phenotype.volume.fluid / 
-			(  phenotype.volume.total + 1e-16 ); 
-		pCell_to_eat->phenotype.volume.fluid_fraction = 0.0; 
+			phenotype.volume.fluid_fraction = phenotype.volume.fluid / 
+				(  phenotype.volume.total + 1e-16 ); 
+			pCell_to_eat->phenotype.volume.fluid_fraction = 0.0; 
 
-		phenotype.volume.cytoplasmic_to_nuclear_ratio = phenotype.volume.cytoplasmic_solid / 
-			( phenotype.volume.nuclear_solid + 1e-16 );
+			phenotype.volume.cytoplasmic_to_nuclear_ratio = phenotype.volume.cytoplasmic_solid / 
+				( phenotype.volume.nuclear_solid + 1e-16 );
+				
+			// update corresponding BioFVM parameters (self-consistency) 
+			set_total_volume( phenotype.volume.total ); 
+			pCell_to_eat->set_total_volume( 0.0 ); 
 			
-		// update corresponding BioFVM parameters (self-consistency) 
-		set_total_volume( phenotype.volume.total ); 
-		pCell_to_eat->set_total_volume( 0.0 ); 
-		
-		// absorb the internalized substrates 
-		
-		// multiply by the fraction that is supposed to be ingested (for each substrate) 
+			// absorb the internalized substrates 
+			
+			// multiply by the fraction that is supposed to be ingested (for each substrate) 
 
-		*(pCell_to_eat->internalized_substrates) *= 
-			*(pCell_to_eat->fraction_transferred_when_ingested); // 
+			*(pCell_to_eat->internalized_substrates) *= 
+				*(pCell_to_eat->fraction_transferred_when_ingested); // 
 
-		*internalized_substrates += *(pCell_to_eat->internalized_substrates); 
-		static int n_substrates = internalized_substrates->size(); 
-		pCell_to_eat->internalized_substrates->assign( n_substrates , 0.0 ); 	
+			*internalized_substrates += *(pCell_to_eat->internalized_substrates); 
+			static int n_substrates = internalized_substrates->size(); 
+			pCell_to_eat->internalized_substrates->assign( n_substrates , 0.0 ); 	
 
-		// conserved quantitites in custom data during phagocytosis
-		// so that phagocyte cell absorbs the full amount from the engulfed cell;
-		for( int nn = 0 ; nn < custom_data.variables.size() ; nn++ )
-		{
-			if( custom_data.variables[nn].conserved_quantity == true )
+			// conserved quantitites in custom data during phagocytosis
+			// so that phagocyte cell absorbs the full amount from the engulfed cell;
+			for( int nn = 0 ; nn < custom_data.variables.size() ; nn++ )
 			{
-				custom_data.variables[nn].value += 
-				pCell_to_eat->custom_data.variables[nn].value; 			
+				if( custom_data.variables[nn].conserved_quantity == true )
+				{
+					custom_data.variables[nn].value += 
+					pCell_to_eat->custom_data.variables[nn].value; 			
+				}
 			}
-		}
-		for( int nn = 0 ; nn < custom_data.vector_variables.size() ; nn++ )
-		{
-			if( custom_data.vector_variables[nn].conserved_quantity == true )
+			for( int nn = 0 ; nn < custom_data.vector_variables.size() ; nn++ )
 			{
-				custom_data.vector_variables[nn].value += 
-				pCell_to_eat->custom_data.vector_variables[nn].value; 
+				if( custom_data.vector_variables[nn].conserved_quantity == true )
+				{
+					custom_data.vector_variables[nn].value += 
+					pCell_to_eat->custom_data.vector_variables[nn].value; 
+				}
 			}
-		}
-		
-		// trigger removal from the simulation 
-		// pCell_to_eat->die(); // I don't think this is safe if it's in an OpenMP loop 
-		
-		// flag it for removal 
-		// pCell_to_eat->flag_for_removal(); 
+			
+			// trigger removal from the simulation 
+			// pCell_to_eat->die(); // I don't think this is safe if it's in an OpenMP loop 
+			
+			// flag it for removal 
+			// pCell_to_eat->flag_for_removal(); 
 
-		// remove all adhesions 
-		// pCell_to_eat->remove_all_attached_cells();
+			// remove all adhesions 
+			// pCell_to_eat->remove_all_attached_cells();
+		}
 		
 	}
 
-	// things that have their own thread safety 
+	// Teardown is deferred rather than done here. flag_for_removal() queues this
+	// cell, and the serial cells_ready_to_die -> die() -> delete_cell() pass does
+	// the same four steps off-thread. Doing them here means mutating other cells'
+	// lists from inside the mechanics parallel-for, where remove_all_* walks a
+	// vector unlocked while attach/detach_cell_as_spring edits it under the
+	// critical -- a real race, not a theoretical one. 
 	pCell_to_eat->flag_for_removal();
-	pCell_to_eat->remove_all_attached_cells();
-	pCell_to_eat->remove_all_spring_attachments();
 	
 	return; 
 }
@@ -1508,140 +1533,148 @@ void Cell::attack_cell( Cell* pCell_to_attack , double dt )
 
 void Cell::fuse_cell( Cell* pCell_to_fuse )
 {
-	// don't ingest a cell that's already fused or fuse self 
-	if( pCell_to_fuse->phenotype.volume.total < 1e-15 || this == pCell_to_fuse )
+	// don't fuse self 
+	if( this == pCell_to_fuse )
 	{ return; } 
-		
+	
 	// make this thread safe 
 	#pragma omp critical
 	{
+		// don't fuse a cell that's already fused or ingested
+		if( pCell_to_fuse->phenotype.volume.total >= 1e-15)
+		{ 
 
-		// set new position at center of volume 
-			// x_new = (vol_B * x_B + vol_S * x_S ) / (vol_B + vol_S )
-		
-		std::vector<double> new_position = position; // x_B
-		new_position *= phenotype.volume.total; // vol_B * x_B 
-		double total_volume = phenotype.volume.total; 
-		total_volume += pCell_to_fuse->phenotype.volume.total ;  
-
-		axpy( &new_position , pCell_to_fuse->phenotype.volume.total , pCell_to_fuse->position ); // vol_B*x_B + vol_S*x_S
-		new_position /= total_volume; // (vol_B*x_B+vol_S*x_S)/(vol_B+vol_S);
-
-		static double xL = get_default_microenvironment()->mesh.bounding_box[0];		 
-		static double xU = get_default_microenvironment()->mesh.bounding_box[3]; 
+			// set new position at center of volume 
+				// x_new = (vol_B * x_B + vol_S * x_S ) / (vol_B + vol_S )
+			
+			std::vector<double> new_position = position; // x_B
+			new_position *= phenotype.volume.total; // vol_B * x_B 
+			double total_volume = phenotype.volume.total; 
+			total_volume += pCell_to_fuse->phenotype.volume.total ;  
 
-		static double yL = get_default_microenvironment()->mesh.bounding_box[1];		 
-		static double yU = get_default_microenvironment()->mesh.bounding_box[4]; 
+			axpy( &new_position , pCell_to_fuse->phenotype.volume.total , pCell_to_fuse->position ); // vol_B*x_B + vol_S*x_S
+			new_position /= total_volume; // (vol_B*x_B+vol_S*x_S)/(vol_B+vol_S);
 
-		static double zL = get_default_microenvironment()->mesh.bounding_box[2];		 
-		static double zU = get_default_microenvironment()->mesh.bounding_box[5]; 
+			static double xL = get_default_microenvironment()->mesh.bounding_box[0];		 
+			static double xU = get_default_microenvironment()->mesh.bounding_box[3]; 
 
-		if( new_position[0] < xL || new_position[0] > xU || 
-		    new_position[1] < yL || new_position[1] > yU || 
-			new_position[2] < zL || new_position[2] > zU )
-		{
-			std::cout << "cell fusion at " << new_position << " violates domain bounds" << std::endl; 
-			std::cout << get_default_microenvironment()->mesh.bounding_box << std::endl << std::endl; 
-		}
-		position = new_position; 
-		update_voxel_in_container();
+			static double yL = get_default_microenvironment()->mesh.bounding_box[1];		 
+			static double yU = get_default_microenvironment()->mesh.bounding_box[4]; 
 
-		// set number of nuclei 
+			static double zL = get_default_microenvironment()->mesh.bounding_box[2];		 
+			static double zU = get_default_microenvironment()->mesh.bounding_box[5]; 
 
-		state.number_of_nuclei += pCell_to_fuse->state.number_of_nuclei; 
+			if( new_position[0] < xL || new_position[0] > xU || 
+				new_position[1] < yL || new_position[1] > yU || 
+				new_position[2] < zL || new_position[2] > zU )
+			{
+				std::cout << "cell fusion at " << new_position << " violates domain bounds" << std::endl; 
+				std::cout << get_default_microenvironment()->mesh.bounding_box << std::endl << std::endl; 
+			}
+			position = new_position; 
+			update_voxel_in_container();
 
-		// absorb all the volume(s)
+			// set number of nuclei 
 
-		// absorb fluid volume (all into the cytoplasm) 
-		phenotype.volume.cytoplasmic_fluid += pCell_to_fuse->phenotype.volume.cytoplasmic_fluid; 
-		pCell_to_fuse->phenotype.volume.cytoplasmic_fluid = 0.0; 
+			state.number_of_nuclei += pCell_to_fuse->state.number_of_nuclei; 
 
-		phenotype.volume.nuclear_fluid += pCell_to_fuse->phenotype.volume.nuclear_fluid; 
-		pCell_to_fuse->phenotype.volume.nuclear_fluid = 0.0; 
+			// absorb all the volume(s)
 
-		// absorb nuclear and cyto solid volume (into the cytoplasm) 
-		phenotype.volume.cytoplasmic_solid += pCell_to_fuse->phenotype.volume.cytoplasmic_solid; 
-		pCell_to_fuse->phenotype.volume.cytoplasmic_solid = 0.0; 
-		
-		phenotype.volume.nuclear_solid += pCell_to_fuse->phenotype.volume.nuclear_solid; 
-		pCell_to_fuse->phenotype.volume.nuclear_solid = 0.0; 
+			// absorb fluid volume (all into the cytoplasm) 
+			phenotype.volume.cytoplasmic_fluid += pCell_to_fuse->phenotype.volume.cytoplasmic_fluid; 
+			pCell_to_fuse->phenotype.volume.cytoplasmic_fluid = 0.0; 
 
-		// consistency calculations 
-		
-		phenotype.volume.fluid = phenotype.volume.nuclear_fluid + 
-			phenotype.volume.cytoplasmic_fluid; 
-		pCell_to_fuse->phenotype.volume.fluid = 0.0; 
-		
-		phenotype.volume.solid = phenotype.volume.cytoplasmic_solid + 
-			phenotype.volume.nuclear_solid; 
-		pCell_to_fuse->phenotype.volume.solid = 0.0; 
-		
-		phenotype.volume.nuclear = phenotype.volume.nuclear_fluid + 
-			phenotype.volume.nuclear_solid; 
-		pCell_to_fuse->phenotype.volume.nuclear = 0.0; 
+			phenotype.volume.nuclear_fluid += pCell_to_fuse->phenotype.volume.nuclear_fluid; 
+			pCell_to_fuse->phenotype.volume.nuclear_fluid = 0.0; 
 
-		phenotype.volume.cytoplasmic = phenotype.volume.cytoplasmic_fluid + 
-			phenotype.volume.cytoplasmic_solid; 
-		pCell_to_fuse->phenotype.volume.cytoplasmic = 0.0; 
-		
-		phenotype.volume.total = phenotype.volume.nuclear + 
-			phenotype.volume.cytoplasmic; 
-		pCell_to_fuse->phenotype.volume.total = 0.0; 
+			// absorb nuclear and cyto solid volume (into the cytoplasm) 
+			phenotype.volume.cytoplasmic_solid += pCell_to_fuse->phenotype.volume.cytoplasmic_solid; 
+			pCell_to_fuse->phenotype.volume.cytoplasmic_solid = 0.0; 
+			
+			phenotype.volume.nuclear_solid += pCell_to_fuse->phenotype.volume.nuclear_solid; 
+			pCell_to_fuse->phenotype.volume.nuclear_solid = 0.0; 
 
-		phenotype.volume.fluid_fraction = phenotype.volume.fluid / 
-			(  phenotype.volume.total + 1e-16 ); 
-		pCell_to_fuse->phenotype.volume.fluid_fraction = 0.0; 
+			// consistency calculations 
+			
+			phenotype.volume.fluid = phenotype.volume.nuclear_fluid + 
+				phenotype.volume.cytoplasmic_fluid; 
+			pCell_to_fuse->phenotype.volume.fluid = 0.0; 
+			
+			phenotype.volume.solid = phenotype.volume.cytoplasmic_solid + 
+				phenotype.volume.nuclear_solid; 
+			pCell_to_fuse->phenotype.volume.solid = 0.0; 
+			
+			phenotype.volume.nuclear = phenotype.volume.nuclear_fluid + 
+				phenotype.volume.nuclear_solid; 
+			pCell_to_fuse->phenotype.volume.nuclear = 0.0; 
 
-		phenotype.volume.cytoplasmic_to_nuclear_ratio = phenotype.volume.cytoplasmic_solid / 
-			( phenotype.volume.nuclear_solid + 1e-16 );
+			phenotype.volume.cytoplasmic = phenotype.volume.cytoplasmic_fluid + 
+				phenotype.volume.cytoplasmic_solid; 
+			pCell_to_fuse->phenotype.volume.cytoplasmic = 0.0; 
 			
-		// update corresponding BioFVM parameters (self-consistency) 
-		set_total_volume( phenotype.volume.total ); 
-		pCell_to_fuse->set_total_volume( 0.0 ); 
+			phenotype.volume.total = phenotype.volume.nuclear + 
+				phenotype.volume.cytoplasmic; 
+			pCell_to_fuse->phenotype.volume.total = 0.0; 
 
-		// absorb the internalized substrates 
-		
-		*internalized_substrates += *(pCell_to_fuse->internalized_substrates); 
-		static int n_substrates = internalized_substrates->size(); 
-		pCell_to_fuse->internalized_substrates->assign( n_substrates , 0.0 ); 	
+			phenotype.volume.fluid_fraction = phenotype.volume.fluid / 
+				(  phenotype.volume.total + 1e-16 ); 
+			pCell_to_fuse->phenotype.volume.fluid_fraction = 0.0; 
 
-		// set target volume(s)
+			phenotype.volume.cytoplasmic_to_nuclear_ratio = phenotype.volume.cytoplasmic_solid / 
+				( phenotype.volume.nuclear_solid + 1e-16 );
+				
+			// update corresponding BioFVM parameters (self-consistency) 
+			set_total_volume( phenotype.volume.total ); 
+			pCell_to_fuse->set_total_volume( 0.0 ); 
 
-		phenotype.volume.target_solid_cytoplasmic += pCell_to_fuse->phenotype.volume.target_solid_cytoplasmic;
-		phenotype.volume.target_solid_nuclear += pCell_to_fuse->phenotype.volume.target_solid_nuclear;
-		
-		// trigger removal from the simulation 
-		// pCell_to_eat->die(); // I don't think this is safe if it's in an OpenMP loop 
-		
-		// flag it for removal 
-		// pCell_to_eat->flag_for_removal(); 
-		// mark it as dead 
-		pCell_to_fuse->phenotype.death.dead = true; 
-		// set secretion and uptake to zero 
-		pCell_to_fuse->phenotype.secretion.set_all_secretion_to_zero( );  
-		pCell_to_fuse->phenotype.secretion.set_all_uptake_to_zero( ); 
-		
-		// deactivate all custom function 
-		pCell_to_fuse->functions.custom_cell_rule = NULL; 
-		pCell_to_fuse->functions.update_phenotype = NULL; 
-		pCell_to_fuse->functions.contact_function = NULL; 
-		pCell_to_fuse->functions.cell_division_function = NULL; 
-		pCell_to_fuse->functions.cell_division_direction_function = NULL; 
-		pCell_to_fuse->functions.volume_update_function = NULL; 
+			// absorb the internalized substrates 
+			
+			*internalized_substrates += *(pCell_to_fuse->internalized_substrates); 
+			static int n_substrates = internalized_substrates->size(); 
+			pCell_to_fuse->internalized_substrates->assign( n_substrates , 0.0 ); 	
 
-		// remove all adhesions 
-		// pCell_to_eat->remove_all_attached_cells();
-		
-		// set cell as unmovable and non-secreting 
-		pCell_to_fuse->is_movable = false; 
-		pCell_to_fuse->is_active = false; 
+			// set target volume(s)
+
+			phenotype.volume.target_solid_cytoplasmic += pCell_to_fuse->phenotype.volume.target_solid_cytoplasmic;
+			phenotype.volume.target_solid_nuclear += pCell_to_fuse->phenotype.volume.target_solid_nuclear;
+			
+			// trigger removal from the simulation 
+			// pCell_to_eat->die(); // I don't think this is safe if it's in an OpenMP loop 
+			
+			// flag it for removal 
+			// pCell_to_eat->flag_for_removal(); 
+			// mark it as dead 
+			pCell_to_fuse->phenotype.death.dead = true; 
+			// set secretion and uptake to zero 
+			pCell_to_fuse->phenotype.secretion.set_all_secretion_to_zero( );  
+			pCell_to_fuse->phenotype.secretion.set_all_uptake_to_zero( ); 
+			
+			// deactivate all custom function 
+			pCell_to_fuse->functions.custom_cell_rule = NULL; 
+			pCell_to_fuse->functions.update_phenotype = NULL; 
+			pCell_to_fuse->functions.contact_function = NULL; 
+			pCell_to_fuse->functions.cell_division_function = NULL; 
+			pCell_to_fuse->functions.cell_division_direction_function = NULL; 
+			pCell_to_fuse->functions.volume_update_function = NULL; 
+
+			// remove all adhesions 
+			// pCell_to_eat->remove_all_attached_cells();
+			
+			// set cell as unmovable and non-secreting 
+			pCell_to_fuse->is_movable = false; 
+			pCell_to_fuse->is_active = false; 
+		}
 
 	}
 
 	// things that have their own thread safety 
+	// Teardown is deferred rather than done here. flag_for_removal() queues this
+	// cell, and the serial cells_ready_to_die -> die() -> delete_cell() pass does
+	// the same four steps off-thread. Doing them here means mutating other cells'
+	// lists from inside the mechanics parallel-for, where remove_all_* walks a
+	// vector unlocked while attach/detach_cell_as_spring edits it under the
+	// critical -- a real race, not a theoretical one. 
 	pCell_to_fuse->flag_for_removal();
-	pCell_to_fuse->remove_all_attached_cells();
-	pCell_to_fuse->remove_all_spring_attachments();
 
 	return; 
 }
@@ -1652,12 +1685,12 @@ void Cell::lyse_cell( void )
 	if( phenotype.volume.total < 1e-15 )
 	{ return; } 	
 	
+	// mark it as dead (before flagging, to match ingest_cell and fuse_cell)
+	phenotype.death.dead = true;
+
 	// flag for removal 
 	flag_for_removal(); // should be safe now 
 	
-	// mark it as dead 
-	phenotype.death.dead = true; 
-	
 	// set secretion and uptake to zero 
 	phenotype.secretion.set_all_secretion_to_zero( );  
 	phenotype.secretion.set_all_uptake_to_zero( ); 
@@ -1671,7 +1704,12 @@ void Cell::lyse_cell( void )
 	
 	// remove all adhesions 
 	
-	remove_all_attached_cells(); 
+	// Teardown is deferred rather than done here. flag_for_removal() queues this
+	// cell, and the serial cells_ready_to_die -> die() -> delete_cell() pass does
+	// the same four steps off-thread. Doing them here means mutating other cells'
+	// lists from inside the mechanics parallel-for, where remove_all_* walks a
+	// vector unlocked while attach/detach_cell_as_spring edits it under the
+	// critical -- a real race, not a theoretical one. 
 	
 	// set volume to zero 
 	set_total_volume( 0.0 ); 
@@ -3457,6 +3495,78 @@ void Cell::remove_all_spring_attachments( void )
 	return; 
 }
 
+void Cell::add_attacker( Cell* pAddMe )
+{
+	#pragma omp critical
+	{
+		if( std::find( state.attacked_by.begin() , state.attacked_by.end() , pAddMe )
+			== state.attacked_by.end() )
+		{ state.attacked_by.push_back( pAddMe ); }
+	}
+	return;
+}
+
+void Cell::remove_attacker( Cell* pRemoveMe )
+{
+	#pragma omp critical
+	{
+		auto search = std::find( state.attacked_by.begin() , state.attacked_by.end() , pRemoveMe );
+		if( search != state.attacked_by.end() )
+		{
+			// copy last entry over it, then shrink by one
+			*search = state.attacked_by.back();
+			state.attacked_by.pop_back();
+		}
+	}
+	return;
+}
+
+void Cell::remove_all_attackers( void )
+{
+	// Swap the list out under the same lock add_attacker / remove_attacker use, then
+	// walk the local copy. Holding the lock across the loop is not an option:
+	// detach_cells_as_spring() takes the same unnamed critical and OpenMP criticals
+	// are not reentrant, so that would deadlock.
+	std::vector<Cell*> attackers;
+	#pragma omp critical
+	{ attackers.swap( state.attacked_by ); }
+
+	for( int i = 0; i < attackers.size() ; i++ )
+	{
+		Cell* pAttacker = attackers[i];
+		if( pAttacker->phenotype.cell_interactions.pAttackTarget == this )
+		{ pAttacker->phenotype.cell_interactions.pAttackTarget = NULL; }
+		detach_cells_as_spring( pAttacker , this );
+	}
+	return;
+}
+
+bool Cell::is_attacked_by( Cell* pCell )
+{
+	bool found = false;
+	#pragma omp critical
+	{
+		found = std::find( state.attacked_by.begin() , state.attacked_by.end() , pCell )
+			!= state.attacked_by.end();
+	}
+	return found;
+}
+
+void Cell::remove_self_from_attacked( void )
+{
+	Cell* pTarget = phenotype.cell_interactions.pAttackTarget;
+	if( pTarget == NULL )
+	{ return; }
+	phenotype.cell_interactions.pAttackTarget = NULL;
+	pTarget->remove_attacker( this );
+	// mutual attackers share one spring, so only drop it once both are done.
+	// "does pTarget attack me?" is asked of our own attacked_by, which we can lock,
+	// rather than by reading pTarget's pAttackTarget unsynchronised.
+	if( !is_attacked_by( pTarget ) )
+	{ detach_cells_as_spring( this , pTarget ); }
+	return;
+}
+
 
 void attach_cells( Cell* pCell_1, Cell* pCell_2 )
 {
@@ -3483,6 +3593,39 @@ void detach_cells_as_spring( Cell* pCell_1 , Cell* pCell_2 )
 {
 	pCell_1->detach_cell_as_spring( pCell_2 );
 	pCell_2->detach_cell_as_spring( pCell_1 );
+	return;
+}
+
+void begin_attack( Cell* pAttacker , Cell* pTarget )
+{
+	// no self-attack, and nothing to do without both ends
+	if( pAttacker == NULL || pTarget == NULL || pAttacker == pTarget )
+	{ return; }
+
+	// a cell attacks at most one target at a time. if it is already attacking
+	// something, close that attack out cleanly rather than orphaning the link.
+	if( pAttacker->phenotype.cell_interactions.pAttackTarget != NULL )
+	{ pAttacker->remove_self_from_attacked(); }
+
+	pAttacker->phenotype.cell_interactions.pAttackTarget = pTarget;
+	pTarget->add_attacker( pAttacker );
+	// spring-link these cells
+	attach_cells_as_spring( pAttacker , pTarget );
+	return;
+}
+
+void end_attack( Cell* pAttacker , Cell* pTarget )
+{
+	if( pAttacker == NULL || pTarget == NULL )
+	{ return; }
+
+	if( pAttacker->phenotype.cell_interactions.pAttackTarget == pTarget )
+	{ pAttacker->phenotype.cell_interactions.pAttackTarget = NULL; }
+	pTarget->remove_attacker( pAttacker );
+	// two cells can attack each other, and attach_cell_as_spring dedupes, so the
+	// two attacks share ONE spring. Only drop it once nobody still needs it.
+	if( !pAttacker->is_attacked_by( pTarget ) )
+	{ detach_cells_as_spring( pAttacker , pTarget ); }
 	return; 
 }
 
diff --git a/core/PhysiCell_cell.h b/core/PhysiCell_cell.h
index 5b7c52d96..4c5e17f6e 100644
--- a/core/PhysiCell_cell.h
+++ b/core/PhysiCell_cell.h
@@ -145,6 +145,10 @@ class Cell_State
 	std::vector<Cell*> attached_cells; 
 	std::vector<Cell*> spring_attachments; 
 
+	// reverse of phenotype.cell_interactions.pAttackTarget. in state, not
+	// phenotype, so it survives convert_to_cell_definition()
+	std::vector<Cell*> attacked_by;
+
 	std::vector<Cell*> neighbors; 
 	std::vector<double> orientation;
 	
@@ -239,6 +243,13 @@ class Cell : public Basic_Agent
 	void detach_cell_as_spring( Cell* pRemoveMe ); // done 
 	void remove_all_spring_attachments( void ); // done 
 
+	void add_attacker( Cell* pAddMe ); // pAddMe has started attacking me
+	void remove_attacker( Cell* pRemoveMe ); // pRemoveMe has stopped attacking me
+	// call these wherever remove_all_spring_attachments() is called
+	void remove_all_attackers( void ); // everyone attacking me stops
+	void remove_self_from_attacked( void ); // I stop attacking whomever I attack
+	bool is_attacked_by( Cell* pCell ); // is pCell currently attacking me?
+
 	// I want to eventually deprecate this, by ensuring that 
 	// critical BioFVM and PhysiCell data elements are synced when they are needed 
 	
@@ -296,6 +307,11 @@ void detach_cells( Cell* pCell_1 , Cell* pCell_2 );
 void attach_cells_as_spring( Cell* pCell_1, Cell* pCell_2 );
 void detach_cells_as_spring( Cell* pCell_1 , Cell* pCell_2 );
 
+// maintain pAttackTarget, attacked_by, and the spring together. use these
+// rather than writing pAttackTarget directly.
+void begin_attack( Cell* pAttacker , Cell* pTarget );
+void end_attack( Cell* pAttacker , Cell* pTarget );
+
 
 std::vector<Cell*> find_nearby_cells( Cell* pCell ); // new in 1.8.0
 std::vector<Cell*> find_nearby_interacting_cells( Cell* pCell ); // new in 1.8.0
diff --git a/core/PhysiCell_standard_models.cpp b/core/PhysiCell_standard_models.cpp
index 4d170a406..2a57d44a7 100644
--- a/core/PhysiCell_standard_models.cpp
+++ b/core/PhysiCell_standard_models.cpp
@@ -1265,7 +1265,9 @@ void standard_cell_cell_interactions( Cell* pCell, Phenotype& phenotype, double
 			{
 				if( UniformRandom() < probability ) 
 				{				
-					pCell->phenotype.cell_interactions.pAttackTarget = pTarget; 
+					// sets pAttackTarget, registers pCell in pTarget's
+					// attacked_by, and spring-links the two
+					begin_attack( pCell , pTarget );
 					attacked = true; 
 					/*					
 					std::cout << "*********   *********  ********  start atack **** " << PhysiCell_globals.current_time << std::endl; 
@@ -1273,8 +1275,6 @@ void standard_cell_cell_interactions( Cell* pCell, Phenotype& phenotype, double
 					<< "attack duration: " << pCell->phenotype.cell_interactions.attack_duration << " "  
 					<< "attack damage rate: " << pCell->phenotype.cell_interactions.attack_damage_rate <<  std::endl; 
 					*/
-					// spring-link these cells 
-					attach_cells_as_spring(pCell,pTarget); 
 				} 
 			}
 
@@ -1316,11 +1316,14 @@ void standard_cell_cell_interactions( Cell* pCell, Phenotype& phenotype, double
 
 	// move effector attack here. 
 
-		if( pCell->phenotype.cell_interactions.pAttackTarget != NULL ) 
+		// read the target once. attack_cell() dereferences without a null check, so
+		// testing this field and then reloading it would let any concurrent clear
+		// hand it a NULL. One load is free and does not depend on which teardown
+		// paths currently run inside this loop. 
+		Cell* pAttackTarget = pCell->phenotype.cell_interactions.pAttackTarget; 
+		if( pAttackTarget != NULL ) 
 		{
-			Cell* pTarget = pCell->phenotype.cell_interactions.pAttackTarget; 
-
-			pCell->attack_cell(pTarget,dt); 
+			pCell->attack_cell(pAttackTarget,dt); 
 			attacked = true; // attacked at least one cell in this time step 
 
 			// attack_cell
@@ -1330,20 +1333,18 @@ void standard_cell_cell_interactions( Cell* pCell, Phenotype& phenotype, double
 			probability = dt / (1e-15 + pCell->phenotype.cell_interactions.attack_duration); 
 
 
-			if( UniformRandom() < probability || pTarget->phenotype.death.dead ) 
+			if( UniformRandom() < probability || pAttackTarget->phenotype.death.dead ) 
 			{
 				/*
 				std::cout << "*********   *********  ********  attack done **** " << PhysiCell_globals.current_time << " " 
 				<< probability << " "
-				<< "attack time: " << pTarget->state.total_attack_time << " " 		
-				<< "damage: " << pTarget->phenotype.cell_integrity.damage <<  " " 		
-				<< "dead? " << (int) pTarget->phenotype.death.dead << " " 
+				<< "attack time: " << pAttackTarget->state.total_attack_time << " " 		
+				<< "damage: " << pAttackTarget->phenotype.cell_integrity.damage <<  " " 		
+				<< "dead? " << (int) pAttackTarget->phenotype.death.dead << " " 
 				<< "damage delivered: " << pCell->phenotype.cell_interactions.total_damage_delivered << std::endl; 
 				*/
 
-				detach_cells_as_spring(pCell,pTarget); 
-
-				pCell->phenotype.cell_interactions.pAttackTarget = NULL; 
+				end_attack( pCell , pAttackTarget );
 			} 
 		} 
 
diff --git a/modules/PhysiCell_MultiCellDS.cpp b/modules/PhysiCell_MultiCellDS.cpp
index 40aac05e9..6ac94df5f 100644
--- a/modules/PhysiCell_MultiCellDS.cpp
+++ b/modules/PhysiCell_MultiCellDS.cpp
@@ -2196,7 +2196,10 @@ int recreate_sim_state(std::string filename, Microenvironment& M,
             {
                 if (cell->ID == pair.second)
                 {
-                    (pair.first)->phenotype.cell_interactions.pAttackTarget = cell;
+                    // restore all three records, not just the pointer: the spring
+                    // is not serialised, so without this a resumed run starts with
+                    // an attack that has no spring behind it
+                    begin_attack( pair.first , cell );
                     if (debug_print)
                     { std::cout << "    cell ID=" << (pair.first)->ID << " attacking  cell ID=" << cell->ID << std::endl; }
                     break;