Repositorio Institucional de Documentos

Trabajo Fin de Máster Máster en Ingeniería Biomédica Cellular response due to substrate stiffness variations: A phenomenological model Autor Aarón Xerach Herrera Martín Directores José Manuel García Aznar Carlos Borau Zamora Departamento de Ingeniería Mecánica Escuela de Ingeniería y Arquitectura Universidad de Zaragoza Septiembre 2013 Cellular response due to substrate stiffness variations: A phenomenological model Appendix I – Calculation report 29 Appendix I - Calculation report 1. Introduction. In this section are exposed the whole results obtained during the development of the work. The main results are included and commented in the Main Report. 2. Checking the model. The first step before start the calculations is checking the model with simplified geometry and boundary conditions. The cubic geometry, composed of eight C3D8 elements, has symmetric like boundary conditions. Therefore, there is no force applied, allowing the geometry to contract freely. Taking into account the expressions developed in the Main Report, the force exerted by the cell can be written as follows: (2.1) knowing that and rearranging the expression, it is obtained the equation below: ( ) (2.2) Since there is no force in this model, the value of can be obtained from the expression above, corresponding its value to = -0.3174. This result is in agreement with the obtained in the computational model showed in figure I.1. Figure I.1. Simplified model employed to check the behaviour equations. The result of the strain is in agreement with the theoretical value: 0.3120 versus 0.3174. Cellular response due to substrate stiffness variations: A phenomenological model. 30 3. Modelling the Mitrossilis et al. experiment. The developed model aims to reproduce the experiment of Mitrossilis et al. using the Finite Element Method. Thus, the model comprises three parts (Figure I.2 and I.3):  Flexible top plate. It has 311.25 μm length, 100 μm width and it is discretized with 1240 S4 elements.  Cell. It has 30 μm height and between 24 μm and 30 μm width depending on the cell shape: Cylinder like or barrel like. It is composed of 4560 C3D8 elements.  Lower plate, modelled as a rigid solid with 100 R3D4 elements. It has enough surface to allow cell contraction and movement. Figure I.2. Developed model to reproduce the experimental study of Mitrossilis et al. with cylinder like cell. Green: Flexible top plate. Orange: Cell. Blue: Rigid lower plate. Figure I.3. Developed model to reproduce the experimental study of Mitrossilis et al. with barrel like cell. Green: Flexible top plate. Orange: Cell. Blue: Rigid lower plate. The cell properties for the active part are set according to the model described in the Main Report, in section 2.1. Material model. The passive components follow the behaviours listed in section 2.3. Behaviour of the passive component of the cytoskeleton. With these settings and varying the elastic modulus of the top plate to achieve different values of stiffness, several cellular responses can be measured. Appendix I – Calculation report. 31 The boundary conditions applied to the model are:  Right end of the flexible top plate totally fixed.  Cell in contact with top and lower plate. The contact can be without roughness, totally attached or with enough friction to allow cell attaching but also movement.  Lower plate totally fixed. The stiffness of the flexible top plate comprises a range from 15 to 1800 nN/μm to address with the stiffness values employed to determine the cell exerted by the cell. Since contact modelling is a challenging task, the cell shape was simplified to a cylinder with the contact interfaces totally attached to both plates. When the convergence was achieved for a wide variety of cantilever stiffnesses, friction coefficients were tested, aiming to a compromise value that allows cell spreading and enough attaching to the substrate. Then, cell shape was changed to an oval shape, increasing the similitude between the Mitrossilis et al. experiment and the computational model. 3.1. Model application with linear elastic behaviour of the passive component. As first approach to the Mitrossilis et al. experiment the cell was modelled with a linear elastic behaviour of the passive component of the model. Two models were tested: Cylinder like cell and barrel like cell. 3.1.1. Cylinder like cell. During the cell contraction, the shape changes as showed in figure I.4, where it can be observed how the cell shrinks axially and changes its shape from a cylinder to a concave form. Figure I.4. Evolution of cell shape along the contraction process with a cylinder like cell. After several calculations with a wide variety of top plate stiffnesses it can be assumed that the axial strain of the cell is higher when the substrate is less stiff. Nevertheless, applying higher stiffnesses, the axial deformation becomes smaller, increasing the transversal strain. This behaviour can be observed in the following pictures, where five representative results are showed. The stiffness of the lower plate for these cases are, in order of appearance: 15 nN/μm, 100 nN/μm, 176 nN/μm, 200 nN/μm and 750 nN/μm. Cellular response due to substrate stiffness variations: A phenomenological model. 32 Figure I.5. Cylinder like cell displacement with linear elastic passive behaviour. Flexible top plate stiffness: 15 nN/μm. Figure I.6. Cylinder like cell displacement with linear elastic passive behaviour. Flexible top plate stiffness: 100 nN/μm. Figure I.7. Cylinder like cell displacement with linear elastic passive behaviour. Flexible top plate stiffness: 176 nN/μm. Appendix I – Calculation report. 33 Figure I.8. Cylinder like cell displacement with linear elastic passive behaviour. Flexible top plate stiffness: 200 nN/μm. Figure I.9. Cylinder like cell displacement with linear elastic passive behaviour. Flexible top plate stiffness: 750 nN/μm. Additionally, the cell without the top and the lower plate is exposed to clarify the explanation given above. It can be seen how the axial displacement becomes smaller as the stiffness of the flexible top plate increases, in contrast to the transversal displacement. Figure I.10. Transversal (left) and axial (right) cell displacement with a flexible top plate stiffness of 15 nN/ μm. Maximum transversal displacement: 4.292 μm. Maximum axial displacement: 11.490 μm. Cellular response due to substrate stiffness variations: A phenomenological model. 34 Figure I.11. Transversal (left) and axial (right) cell displacement with a flexible top plate stiffness of 176 nN/ μm. Maximum transversal displacement: 5.347 μm. Maximum axial displacement: 6.979 μm. Figure I.12. Transversal (left) and axial (right) cell displacement with a flexible top plate stiffness of 750 nN/ μm. Maximum transversal displacement: 7.420 μm. Maximum axial displacement: 3.870 μm. As expected, the axial displacement of the cell decreases as the stiffness of the flexible top plate increases. However, it can be observed to higher stiffnesses that the axial displacement decreases to such degree that it is smaller than the transversal displacement. This behaviour can be explained in a way that the cell feels itself stiffness and react in consequence to its own solicitations. To clarify this explanation, the figure I.13 is given, where it is shown how the axial displacement (blue) becomes smaller until it reach a maximum. On the other hand, the transversal contraction (red) continues increasing. Appendix I – Calculation report. 35 Figure I.13. Maximum displacements obtained in the cell. Passive component characterized as linear elastic. Cylinder like cell. Blue: Axial displacement. Red: Transversal displacement. Finally, after the development of the calculations with several stiffnesses of the flexible top plate with the cylinder like cell, the forces exerted by the cell are shown in figure I.14. Figure I.14. Force obtained for different stiffness values. Passive component characterized as linear elastic. Cylinder like cell. In spite of the appropriate cell behaviour observed when it changes its shape, a plateau force is not achieved. Although it can be observed a change in the trend of the force exerted by the cell when it reaches 200 nN/μm, it is not enough to stabilize the force to a same value. Therefore, an improvement in the model is needed. The first step is changing to another suitable cell shape. 3 4 5 6 7 8 9 10 11 12 13 4 504 1004 1504 2004 Displacement (μm) Stiffness (nN/μm) Maximum displacements in axial and transversal directions. 100 600 1100 1600 2100 2600 3100 3600 4100 10 210 410 610 810 1010 1210 1410 1610 1810 2010 Force (nN) Stiffness (nN/μm) Linear elastic behaviour of the passive component. Cylinder like cell. Cellular response due to substrate stiffness variations: A phenomenological model. 42 Figure I.31. Transversal (left) and axial (right) cell displacement with a flexible top plate stiffness of 750 nN/ μm. Maximum transversal displacement: 3.892 μm. Maximum axial displacement: 3.663 μm. Similarly to the previous cases, the axial displacement of the cell decreases as the stiffness of the flexible top plate increases. Nevertheless, the transversal displacement tends to stabilize after it reaches the maximum, possibly due to the poroelastic behaviour of the passive side. Although the axial displacement does not follow the same trend, it must be taken into account not only the stiffness of the top plate but also the shape of the cell on the plates that makes the cell has more material in one direction than another. Figure I.32. Maximum displacements obtained in the cell. Passive component characterized as poroelastic. Cylinder like cell. Blue: Axial displacement. Red: Transversal displacement. Finally, after the calculation of several cell models with different top plate stiffnesses, it can be obtained the results showed in figure I.33, where it can be observed the trend of the cell force as the stiffness increases. 1 2 3 4 5 6 7 8 9 10 11 4 504 1004 1504 2004 Displacement (μm) Stiffness (nN/μm) Maximum displacements in axial and transversal directions. Appendix I – Calculation report. 43 Figure I.33. Force obtained for different stiffness values. Passive component characterized as poroelastic. Cylinder like cell. Each point corresponds to each result of the model. The main difference with the models above, in which a linear elastic behaviour of the passive side is employed, is that the cell reaches a plateau force established between 2200 nN and 2400 nN once the stiffness of the flexible top plate is set to 750 nN/μm. Nevertheless, although the force exerted by the cell behaves properly, reaching a maximum force after a specific stiffness, the shape development does not agree with the observed in the experiment, in which the cell shrinks axially and spreads transversally. A feasible explanation of this behaviour is given in the conclusion section of the Main Report. 4. Model application to real cell geometry. On the other hand, taking advantage of the relationship between this University and the CIMA, there was studied a real cell model. The cellular model is based on a breast cancer cell line MDA-MB-231 cultured using microfluidic devices. The cell was imaged under confocal microscopy slice by slice over the Z axis, with 1 μm distance between them. After the segmentation to smooth the cell boundary and to give the cytoskeleton staining, the segmented cell is transformed into a cloud of points through the sampled points obtained at the cell boundary. Then, the cell is discretized into three-dimensional volumes via tetrahedralization, which creates the elements to work with the finite element model. Hence, the cell is characterized in the same way as the previous case, with a passive and an active side, duplicating the elements. 100 600 1100 1600 2100 2600 10 210 410 610 810 1010 1210 1410 1610 1810 2010 Force (nN) Stiffness (nN/μm) Poroelastic behaviour of the passive component. Cylinder like cell. Cellular response due to substrate stiffness variations: A phenomenological model. 44 Figure I.34. Cellular model obtained in the Centro de Investigación Médica Aplicada of Navarra. 4.1. Cell embedded in a substrate. As first approach, the cell is embedded in a cubic substrate of 194 μm side characterized both as linear elastic and hyperelastic behaviour. The values are 1 KPa and a Poisson ratio of 0.49 for the linear elastic model. Then, for the hyperelastic behaviour, the values for the characterization are obtained via uniaxial test. The parameters are given in the table I.1. Figure I.35. Cellular model (green) inside the substrate (blue). Appendix I – Calculation report. 45 σn λ 0.304905 9.9995e-05 0.436859 0.000145989 0.631863 0.000213977 0.913933 0.000312951 1.32549 0.000457895 1.90341 0.000668776 2.76873 0.000977522 4.01796 0.00142898 5.81187 0.00208782 8.43461 0.00305533 12.197 0.00446004 17.5809 0.00651871 25.2385 0.0095245 35.9336 0.0139029 50.3498 0.0202927 68.0394 0.0294617 86.6822 0.0428679 102.604 0.0620354 112.008 0.089475 113.088 0.128393 106.74 0.182322 Table I.1. Hyperelastic parameters, where σn is the nominal stress in kPa and λ is the nominal strain. On the other hand, the cell properties are set as linear elastic with an elastic modulus of 75 KPa and near-incompressible behaviour with a Poisson ratio of 0.49. It must be taken into account that, unlike a real cell, this cell is totally attached to the substrate and, hence, it can exert force around its entire surface. This model has 1074096 C3D4 elements divided as follows:  529648 elements corresponding to the substrate.  272223 elements corresponding to the passive side of the cell.  272223 elements corresponding to the active side of the cell. Finally, the boundary conditions are applied to the substrate surface, making it totally fixed in order to allow the cell to feel the mechanical properties of its surrounding matrix. Cellular response due to substrate stiffness variations: A phenomenological model. 46 4.1.1. Linear elastic substrate application. After the model application to the cell geometry with the linear elastic substrate, the results are shown in the figures I.36 and I.37. Figure I.36. Side view of the cell contraction obtained after the model application. Undeformed shape (translucent) and deformed shape (inside). Displacement displayed. It can be observed a cell contraction between 0.66 and 0.03 μm. Scale factor: 15. Figure I.37. Side view of the cell contraction obtained after the model application. Undeformed shape (translucent) and deformed shape (inside). Displacement displayed. It can be observed a cell contraction between 0.66 and 0.03 μm. Scale factor: 15. As can be observed in figures I.36 and I.37, cell model involves a whole cell contraction without shape change. Appendix I – Calculation report. 47 4.1.2. Hyperelastic substrate application. After the model application to the cell geometry with the hyperelastic substrate, the results are the exposed in the figures I.38 and I.39. Figure I.38. Side view of the cell contraction obtained after the model application. Undeformed shape (translucent) and deformed shape (inside). Displacement displayed. It can be observed a cell contraction between 0.50 and 0.02 μm, with a protrusion accentuation. Scale factor: 15. The substrate is not displayed. Figure I.39. Rear view of the cell contraction obtained after the model application. Undeformed shape (translucent) and deformed shape (inside). Displacement displayed. It can be observed a cell contraction between 0.50 and 0.02 μm, with a protrusion accentuation. Scale factor: 15. The substrate is not displayed. As shown in the figures I.38 and I.39, cell model involves not only a whole cell contraction but also a protrusion accentuation, in agreement with the probing role of the lamellipodia and filopodia in mechanosensing mechanism. Cellular response due to substrate stiffness variations: A phenomenological model. 48 4.1.3. Comparison between elastic and hyperelastic substrate. A brief comparison of the influence in cell behaviour depending on the properties of the substrate can be observed in the figure I.40. Figure I.40. Rear view of the cell contraction obtained after the model application. Undeformed shape (translucent) and deformed shape (inside). Displacement displayed. It can be observed a cell contraction between 0.50 and 0.02 μm, with a protrusion accentuation. Scale factor: 15. The substrate is not displayed. As previously exposed, cell strain differs depending on whether the cell is in a linear elastic substrate or in a hyperelastic substrate, producing not only a whole cell contraction but also an accentuation of the protrusion in the second case. 4.2. Simulating focal adhesions. Since the cells are not totally attached to its substrate, the case developed above must be refined. Thus, the substrate is eliminated and it is replaced by boundary conditions applied to the main protrusion that can be found on the cell surface. These boundary conditions are modelled like springs that simulate the focal adhesion that could make the cell when is embedded in its ECM. In addition, applying several stiffnesses to the springs, different kinds of the substrates can be modelled. In this way, six adhesions around the cell surface are created, employing springs on each nodal surface of the protrusions. The spring stiffnesses range from 1 nN/μm to 2500 nN/μm aiming to simulate different substrates. The position of these adhesions is shown below: Appendix I – Calculation report. 49 Figure I.41. Cell (green) with focal adhesions placed on the protrusions (purple). Taking into account this kind of boundary conditions, the model is composed with a total of 545642 elements divided as follows:  272224 C3D4 elements corresponding to the passive side of the cell.  272224 C3D4 elements corresponding to the active side of the cell.  1194 spring elements corresponding to the focal adhesions. 4.2.1. Cell without nucleus. Once the boundary conditions are set, the cell is modelled, as seen previously, with a linear elastic passive side with an elastic modulus of 1 kPa and a Poisson ratio of 0.45. In the same way, the active side is characterized as linear elastic with an elastic modulus of 10 kPa and a Poisson ratio of 0.45. Following, there are exposed the principal results of the cell model, where can be observed different displacements and stresses depending on the stiffness of the focal adhesion. As expected, the minimum cell displacement occurs in the focal adhesions where it is attached whereas the cell body shrinks while it exerts force. Moreover, the stress distribution is uniform throughout the cell body, increasing in the cell adhesions. Cellular response due to substrate stiffness variations: A phenomenological model. 50 Figure I.42. Displacements obtained with 7.5 nN/μm of focal adhesion stiffness. Figure I.43. Stresses obtained with 7.5 nN/μm of focal adhesion stiffness. Figure I.44. Displacements obtained with 100 nN/μm of focal adhesion stiffness. Appendix I – Calculation report. 51 Figure I.45. Stresses obtained with 100 nN/μm of focal adhesion stiffness. Figure I.46. Displacements obtained with 1000 nN/μm of focal adhesion stiffness. Figure I.47. Stresses obtained with 1000 nN/μm of focal adhesion stiffness. Cellular response due to substrate stiffness variations: A phenomenological model. 58 Figure I.62. Stresses obtained with 100 nN/μm of focal adhesion stiffness. Figure I.63. Displacements obtained with 1000 nN/μm of focal adhesion stiffness. Figure I.64. Stresses obtained with 1000 nN/μm of focal adhesion stiffness. The main difference that can be found between both models with nucleus is the solicitation of this part. In the second case, the nucleus is stiffer and hence, dissimilarities can be found in the stress field. Appendix I – Calculation report. 59 Attending to focal adhesions and cell displacements, there is little difference between them. However, as can be seen in figure I.65, less displacement is achieved by this later model. Figure I.65. Focal adhesion displacement (red) and cell displacement (blue) depending on the stiffness of the substrate. It is remarkable that the displacement of the focal adhesions and the cell with lower stiffnesses are greater than the cases above. This behaviour can be explained because the nucleus is stiffer than the previous case and not only feels the stiffness of the focal adhesions but also the stiffness of its nucleus. In this way, the cell exerts more force and increases its displacements. Nevertheless, the stress generated in the cell is similar to the previous cases as can be seen in figure I.66 The stress reaches a maximum value around 14.5 kPa when the stiffness of the focal adhesion reaches 500 nN/ μm. Figure I.66. Maximum stress obtained within the cell with nucleus depending on the stiffness of the simulated adhesions. The elastic modulus of the nucleus is set to 100 kPa. 5.64 5.66 5.68 5.7 5.72 5.74 5.76 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 110 100 1000 10000 Blue: Cell displacement (μm) Red: Focal adhesion dispalcement (μm) Focal adhesion stiffness (nN/μm) Focal adhesion and cell displacement 6000 7000 8000 9000 10000 11000 12000 13000 14000 15000 0 500 1000 1500 2000 2500 Maximum stress (Pa) Focal adhesion stiffness (nN/μm) Maximum stress in the cell with a nucleus elastic modulus of 100 kPa Cellular response due to substrate stiffness variations: A phenomenological model Appendix II – Programming codes 63 Appendix II – Programing codes 1. Introduction. In this appendix is shown the main programming codes employed for the model development that comprises three codes: The phenomenological model, the element duplication subroutine and the substrate deletion and nucleus creation. 2. The phenomenological model. The phenomenological model is programmed as combination of user subroutines in Abaqus® and then compiled with Intel® Fortran Compiler. The code, with comments, is exposed below. C Variables to be defined. MODULE GLOBVAR PARAMETER (NELEMS=6243) PARAMETER (NNODES=2736) PARAMETER (NPGAUSS=8) REAL*8 EMIN(NELEMS,NPGAUSS) REAL*8 EAVE(NELEMS,NPGAUSS) REAL*8 EMAX(NELEMS,NPGAUSS) REAL*8 EXPANS_ELE(NELEMS,3) REAL*8 EXPANS_LAST(NELEMS,3) REAL*8 EcompMAX REAL*8 EcompAVE REAL*8 EcompMIN END MODULE C -------------------------UEXTERNALDB-------------------------------C C It allows to work with several subroutines. C C--------------------------------------------------------------------C SUBROUTINE UEXTERNALDB(LOP,LRESTART,TIME,DTIME,KSTEP,KINC) USE GLOBVAR INCLUDE 'ABA_PARAM.INC' DIMENSION TIME(2) CHARACTER(256) FILENAME CHARACTER(256) JOBDIR if (LOP.EQ.0) then C LOP 0: the subroutine is called at the start of the analysis. C Values are set to 0. EMIN=0.0d0 EAVE=0.0d0 EMAX=0.0d0 EXPANS_ELE=0.0d0 EXPANS_LAST=0.0d0 endif if (LOP.EQ.1) then Cellular response due to substrate stiffness variations: A phenomenological model. 64 C LOP 1: the subroutine is called at the start of the increment. call GETOUTDIR(JOBDIR,LENJOBDIR) endif if (LOP.EQ.2) then C LOP 2: the subroutine is called at the end of the increment. EXPANS_ELE=EXPANS_LAST endif if (LOP.EQ.3) then C LOP 3: the subroutine iss called at the end of the analysis. close(101) endif RETURN END C -----------------------------URDFIL--------------------------------C C It allows to read the result files. C C--------------------------------------------------------------------C SUBROUTINE URDFIL(LSTOP,LOVRWRT,KSTEP,KINC,DTIME,TIME) USE GLOBVAR INCLUDE 'ABA_PARAM.INC' DIMENSION ARRAY(6000),JRRAY(NPRECD,6000),TIME(2) EQUIVALENCE(ARRAY(1),JRRAY(1,1)) integer KELEM,KPOINT CALL POSFIL(KSTEP,KINC,ARRAY,JRCD) DO K1=1,999999999 CALL DBFILE(0,ARRAY,JRCD) IF (JRCD.NE.0) GO TO 110 KEY=JRRAY(1,2) IF (KEY.EQ.1) THEN C KELEM, element number (Record key 1, record type 1 y 2 + 2, to define JRRAY(1,1). C KPOINT, integration point. KELEM=JRRAY(1,3) KPOINT=JRRAY(1,4) ENDIF C Record KEY 21: Total strain. First, second, etc. C Record KEY 403: Principal strains en orden: MIN, MID, MAX. IF (KEY.EQ.21) THEN EMIN(KELEM,KPOINT)=ARRAY(3) EAVE(KELEM,KPOINT)=ARRAY(4) EMAX(KELEM,KPOINT)=ARRAY(5) write(101,*),EMAX(KELEM,KPOINT),EAVE(KELEM,KPOINT), Appendix II – Programming codes. 65 * EMIN(KELEM,KPOINT) ENDIF ENDDO 110 CONTINUE return end C-----------------------------UEXPAN---------------------------------C C It calculates the element expansion. C C--------------------------------------------------------------------C SUBROUTINE UEXPAN(EXPAN,DEXPANDT,TEMP,TIME,DTIME,PREDEF, 1 DPRED,STATEV,CMNAME,NSTATV,NOEL) USE GLOBVAR INCLUDE 'ABA_PARAM.INC' CHARACTER*80 CMNAME REAL*8 KACT, FMAX, KPAS, INFLIM, SUPLIM, EMAXX, ALPHA(3) REAL*8 TOLERANCIA DIMENSION EXPAN(*),DEXPANDT(*),TEMP(2),TIME(2),PREDEF(*), 1 DPRED(*),STATEV(NSTATV) C KACT: Actin stiffness. C KPAS: Passive stiffness. C FMAX: Maximum force of the contractil system. C INFLIM: Upper strain limit. C SUPLIM: Lower strain limit. KACT=10000.0d0 FMAX=2500.0d0 INFLIM=-0.4d0 SUPLIM=0.4d0 LIM=FMAX/KACT C Average of the Gauus point in each element for EMAX, EAVE, EMIN. EcompMAX=SUM(EMAX(NOEL,:))/NPGAUSS EcompAVE=SUM(EAVE(NOEL,:))/NPGAUSS EcompMIN=SUM(EMIN(NOEL,:))/NPGAUSS C EXPAN(1): X, EXPAN(2): Y, EXPAN(3): Z C------------------EMAX ------------------C IF (EcompMAX .LT. INFLIM)THEN C The strain of the last step is read, that is the total strain. C The new expansion increment will be ALPHA minus the last strain state. ALPHA(2)= EcompMAX EXPAN(2)= ALPHA(2)- EXPANS_ELE(NOEL,2) ELSE IF ((EcompMAX .GE. INFLIM) .AND. (EcompMAX .LE. LIM))THEN C Alpha is calculated (total strain) as function of the last strain state. Cellular response due to substrate stiffness variations: A phenomenological model. 66 ALPHA(2)= ((EcompMAX-(FMAX/KACT))*INFLIM*KACT)/(KACT*INFLIMFMAX) C The new value for the strain increment is de previously calculated minus the element strain value of the previous state. EXPAN(2)= ALPHA(2)- EXPANS_ELE(NOEL,2) ELSE IF ((EcompMAX .GT. LIM) .AND. (EcompMAX .LE. SUPLIM)) THEN ALPHA(2)= ((EcompMAX-(FMAX/KACT))*SUPLIM*KACT)/(KACT*SUPLIMFMAX) EXPAN(2)= ALPHA(2)- EXPANS_ELE(NOEL,2) ELSE IF (EcompMAX .GT. SUPLIM)THEN ALPHA(2)= EcompMAX EXPAN(2)= ALPHA(2)- EXPANS_ELE(NOEL,2) END IF C------------------EAVE------------------C IF (EcompAVE .LT. INFLIM)THEN ALPHA(1)= EcompAVE EXPAN(1)= ALPHA(1)- EXPANS_ELE(NOEL,1) ELSE IF ((EcompAVE .GE. INFLIM) .AND. (EcompAVE .LE. LIM))THEN ALPHA(1)= ((EcompAVE-(FMAX/KACT))*INFLIM*KACT)/(KACT*INFLIMFMAX) EXPAN(1)= ALPHA(1)- EXPANS_ELE(NOEL,1) ELSE IF ((EcompAVE .GT. LIM) .AND. (EcompAVE .LE. SUPLIM)) THEN ALPHA(1)= ((EcompAVE-(FMAX/KACT))*SUPLIM*KACT)/(KACT*SUPLIMFMAX) EXPAN(1)= ALPHA(1)- EXPANS_ELE(NOEL,1) ELSE IF (EcompAVE .GT. SUPLIM)THEN ALPHA(1)= EcompAVE EXPAN(1)= ALPHA(1)- EXPANS_ELE(NOEL,1) END IF C------------------EMIN------------------C IF (EcompMIN .LT. INFLIM)THEN ALPHA(3)= EcompMIN EXPAN(3)= ALPHA(3)- EXPANS_ELE(NOEL,3) ELSE IF ((EcompMIN .GE. INFLIM) .AND. (EcompMIN .LE. LIM))THEN ALPHA(3)= ((EcompMIN-(FMAX/KACT))*INFLIM*KACT)/(KACT*INFLIMFMAX) EXPAN(3)= ALPHA(3)- EXPANS_ELE(NOEL,3) ELSE IF ((EcompMIN .GT. LIM) .AND. (EcompMIN .LE. SUPLIM)) THEN Appendix II – Programming codes. 67 ALPHA(3)= ((EcompMIN-(FMAX/KACT))*SUPLIM*KACT)/(KACT*SUPLIMFMAX) EXPAN(3)= ALPHA(3)- EXPANS_ELE(NOEL,3) ELSE IF (EcompMIN .GT. SUPLIM)THEN ALPHA(3)= EcompMIN EXPAN(3)= ALPHA(3)- EXPANS_ELE(NOEL,3) END IF EXPANS_LAST(NOEL,1)=ALPHA(1) EXPANS_LAST(NOEL,2)=ALPHA(2) EXPANS_LAST(NOEL,3)=ALPHA(3) RETURN END 3. Element duplication. As previously exposed, the model is based on a set of duplicated elements. Nevertheless, depending on the design, the number of elements to be duplicated can be unapproachable, so a Matlab® subroutine is created to facilitate the addition of the duplicated elements. The program is exposed below. clear all load Elset.txt; [Elem, Nodes] = size(Elset); % Elem (rows): Set element number. % Nodos (columns): Node number +1, corresponding to element number. for i=1:1:Elem for j=1:1:Nodes if j==1; Add(i,j)=Elem; else Add(i,j)=0; end end end % First column stores the last element number of the set. % Rest of the columns correspond to the node position of that element % and are set to 0 to avoid modifications in the element shape. % Generated matrix is added to Elset matrix. This way, the elements % are duplicated avoiding node modification. NewElset=Elset+Add; fid=fopen('Elem_Hex.txt','w'); fprintf(fid,'%2d, %2d, %2d, %2d, %2d, %2d, %2d, %2d, %2d \n', NewElset'); fclose(fid); Cellular response due to substrate stiffness variations: A phenomenological model. 68 4. Substrate deletion and nucleus creation. To make the deletion of the substrate is needed a program to eliminate all elements that are not included in the cell. Additionally, to create the nucleus, the cell centroid must be calculated and then, taking it as reference, it is described a sphere of a given radius. All elements within the sphere will form the nucleus while the elements outside will make up the cytoskeleton. clear all clc % Load files nodes_raw=load('Nodos.txt'); conec_raw=load('Conectividad.txt'); nelem=size(conec_raw,1); conec=[conec_raw(:,2) conec_raw(:,3) conec_raw(:,4) conec_raw(:,5)]; node_list=[conec_raw(:,2); conec_raw(:,3); conec_raw(:,4); conec_raw(:,5)]; node_list=sort(node_list); %Sort the list node_list=unique(node_list); %Delete duplicated nodes nodes_coord=nodes_raw(node_list,:); nnodes=size(nodes_coord,1); %Final number of nodes conec_final=conec; c=0; for n=1:nnodes c=c+1; temp_index=nodes_coord(n,1); %the index from the orginal mesh nodes_coord(n,1)=c; %the new index k=find(conec==temp_index);%Find all the nodes numbered as temp_index conec_final(k)=c; %Change them to the new index if mod(c,1000)==0 fprintf('%d nodes from %d processed (~%2.2f%%)... \n',c,nnodes,c/nnodes*100); end end fprintf('--------------------- \n'); fprintf('Printing new files... \n'); fi=fopen('New_coords.txt','w'); fprintf(fi,'%d , %3.6f , %3.6f , %3.6f \n',nodes_coord'); fclose(fi); conec_final=[(1:nelem)' conec_final]; fi=fopen('New_conec.txt','w'); Appendix II – Programming codes. 69 fprintf(fi,'%6d , %6d , %6d , %6d, %6d \n',conec_final'); fclose(fi); fprintf('--------------------- \n'); fprintf('Calculating nucleus... \n'); elem_centroids=zeros(nelem,3); temp_cent=[0 0 0]; % Calculate centroid of each element for n=1:nelem for k=1:4 k_node=conec_final(n,k+1); temp_cent=temp_cent+[nodes_coord(k_node,2) nodes_coord(k_node,3) nodes_coord(k_node,4)]; end elem_centroids(n,:)=temp_cent/4; %Centroid of the element n temp_cent=[0 0 0]; end % Calculate cell centroid cell_centroid=sum(elem_centroids)/nelem ; % Calculate list of elements within a sphere of a specific radius (the nucleus) list_nucleus=zeros(nelem,1); radius=5 for n=1:nelem dist=norm(cell_centroid-elem_centroids(n,:)); if dist<radius list_nucleus(n,1)=n; end if dist>=radius list_cyto(n,1)=n; end end list_nucleus=list_nucleus(list_nucleus~=0); %Clear the zeros fprintf('--------------------- \n'); fprintf('Printing nucleus... \n'); fi=fopen('Nucleus.txt','w'); fprintf(fi,'%5d , %5d, %5d, %5d, %5d, %5d, %5d, %5d, \n',list_nucleus); fclose(fi); fprintf('--------------------- \n'); fprintf('Printing cytoskeleton... \n'); fi=fopen('Cytoskeleton.txt','w');