From 5e2e5f3735894e10f71f0f6e925c08552f9f83b4 Mon Sep 17 00:00:00 2001
From: Jaro Camphuijsen <jjecamphuijsen@gmail.com>
Date: Thu, 4 Mar 2021 12:57:47 +0100
Subject: [PATCH] add more realistic adhesion force options

---
 model_define.h                |  2 +-
 model_routine_mech_intrct.cpp | 62 +++++++++++++++++++++++++++++++----
 2 files changed, 57 insertions(+), 7 deletions(-)

diff --git a/model_define.h b/model_define.h
index f2bcec9..60f1840 100644
--- a/model_define.h
+++ b/model_define.h
@@ -136,7 +136,7 @@ const REAL A_MCARRIER_DENSITY_PER_UB =  1000 / (32.0 *32.0 * 32.0 ) ;  // 2000 /
 const REAL INIT_CELLS_PER_MICROCARRIER = 4; //10;
 const REAL A_CELL_D_MAX[ NUM_AGENT_TYPES ] = { 150.0 * 1.25, 20.0 * 1.25 };
 
-
+const S32 ADHESION_TYPE = 2; //change the type of adhesion implementation, 1 for tanh based, 2 for piecewise linear, 3 for lennart-jones
 const REAL ADHESION_S = 0.01;
 const REAL RANDOM_VIBRATION_SCALE = 0.05;
 
diff --git a/model_routine_mech_intrct.cpp b/model_routine_mech_intrct.cpp
index cb66fd4..732649b 100644
--- a/model_routine_mech_intrct.cpp
+++ b/model_routine_mech_intrct.cpp
@@ -21,6 +21,9 @@ NOTICE:  These data were produced by Battelle Memorial Institute (BATTELLE) unde
 using namespace std;
 
 #if HAS_SPAGENT
+
+static REAL addAdhesionForce( REAL& Fij, const REAL& xij, const REAL& sij);
+
 void ModelRoutine::initJunctionSpAgent( const VIdx& vIdx0, const SpAgent& spAgent0, const UBEnv& ubEnv0, const VIdx& vIdx1, const SpAgent& spAgent1, const UBEnv& ubEnv1, const VReal& vDir/* unit direction vector from spAgent1 to spAgent0 */, const REAL dist, BOOL& link, JunctionEnd& end0/* dummy if link == false */, JunctionEnd& end1/* dummy if link == false */ ) {
 	/* MODEL START */
     link = false;
@@ -72,6 +75,7 @@ void ModelRoutine::computeMechIntrctSpAgent( const S32 iter, const VIdx& vIdx0,
     REAL R1 = A_AGENT_SHOVING_SCALE[type1]*spAgent1.state.getModelReal( CELL_MODEL_REAL_RADIUS ) ;
     REAL dist_threshold = R0 + R1;
     REAL D = R0 + R1 - 0.5*A_AGENT_SHOVING_LIMIT[type0] - 0.5*A_AGENT_SHOVING_LIMIT[type1];
+    REAL Fij = 0.0;
     REAL mag = 0.0;
     REAL stress = 0.0 ; 
     REAL xij  = D - dist ;
@@ -86,7 +90,7 @@ void ModelRoutine::computeMechIntrctSpAgent( const S32 iter, const VIdx& vIdx0,
             }
             else{
                 // compute elastic force
-                REAL Fij = 0.5 * xij * tanh(FABS(xij)*sij);
+                Fij += addAdhesionForce(Fij, xij, sij);
                 mag = mag + Fij ;
                 stress = stress + dist * Fij ;
 
@@ -118,11 +122,13 @@ void ModelRoutine::computeMechIntrctSpAgent( const S32 iter, const VIdx& vIdx0,
                }
 
                // add force rigth away
-               REAL D = R0 + R1;
-               REAL xij  = D - dist  ;
-               REAL Fij = 0.5 * xij * tanh(FABS(xij)*sij);
-               mag = mag + Fij ;
-               stress = stress + dist * Fij;
+                REAL D = R0 + R1;
+                REAL xij  = D - dist  ;
+                REAL Fij = 0.0;
+                Fij += addAdhesionForce(Fij, xij, sij);
+                
+                mag = mag + Fij ;
+                stress = stress + dist * Fij;
             }
         }
         
@@ -150,5 +156,49 @@ void ModelRoutine::computeMechIntrctSpAgent( const S32 iter, const VIdx& vIdx0,
     return;
 
 }
+static REAL addAdhesionForce( REAL& Fij, const REAL& xij, const REAL& sij){
+    // compute elastic force
+    switch (ADHESION_TYPE) {
+
+        //original tanh based adhesion force (more info ask Boris Aguilar)
+        case 1:
+            Fij += 0.5 * xij * tanh(FABS(xij)*sij);
+            break;
+
+        // piecewise linear based on dx.doi.org/10.1021/la500045q
+        // bond strenght sij should be 1 for cell-microcarrier bond in above literature
+        case 2: 
+            if ( xij < 1.5 ) { //distance in micron
+                Fij += sij * xij * 30; //force in nN
+            }
+            else if ( xij < 5) {
+                Fij += sij * 45 - 0.3 * (xij - 1.5);
+            }
+            else if ( xij < 13) {
+                Fij += sij * 43 - 5 * (xij - 5);
+            }
+            else if ( xij < 16.5) {
+                Fij += sij * 3.5 - 0.9 * (xij - 13);
+            }
+            break;
+
+        // adhesion force with Lennard-Jones potential
+        case 3:{
+            REAL sig; 
+            REAL eps = 45 ;
+            REAL LJBase = sig / FABS(xij);
+
+            Fij += 4 * eps * ( pow( LJBase , 6) - pow( LJBase , 12) ) * sij;
+            break;
+        }
+        // adhesion force by integrin simulation
+        case 4:
+            //TODO
+            break;
+
+    }
+
+    return Fij;
+}
 #endif