in the united states , 3 million cases of suspected child abuse are filed each year , but only 30% of these lead to a conviction , indicating that diagnosing child abuse is difficult . the presence of multiple bruises of different color suggests long term abuse , as different colors imply different stages of healing , and thus different times of origin for each bruise . therefore , accurate determination of the age of bruises has the potential to contribute to an improved diagnosis of physical abuse , as it can confirm or exclude the presence of a possible perpetrator at the determined time point . the color of the bruise changes over time [ 18 , 39 ] , hence color has become the choice for age determination . these colors and their corresponding age are documented in standardized protocols available to the physician . however , there is an enormous variability in color descriptions versus age of the bruise which seriously complicates the age determination of the bruise by visual inspection alone [ 3 , 18 , 22 , 26 , 32 , 33 , 39 ] . therefore , it has been attempted to objectively determine the color of a bruise by applying reflectance spectroscopy [ 6 , 14 , 24 , 2830 ] . however , and perhaps more importantly , color is not uniquely characterizing for the age of a bruise since bonhert et al .   showed a relation between color and depth of the bruise : i.e. a deeper bruise has a more blue appearance than a superficial bruise of equal age . in other words , two bruises of equal age but located in different depths in the skin have different appearances . unfortunately , therefore , accurate age determination is currently not possible . following randeberg et al . , we hypothesize that an increased accuracy of age determination of bruises is only possible when the pathophysiology of bruise formation and healing is better understood than now . randeberg et al . constructed an analytic model of formation and healing of bruises , based on darcy s law of convection and fick s law of diffusion , to describe fluid transport in tissue and distribution of blood and hemoglobin breakdown products , and computing the temporal distribution of blood and hemoglobin breakdown products integrated over the whole bruise area . by applying an integrating sphere to illuminate a bruise and measuring the spectrum of the remitted light , randeberg et al the optical data were interpreted with an analytic model and an accuracy of 3  h for a fresh bruise ( less than 1  day old ) , and 1  day for a bruise up to 10  days old was reported . in randeberg s model , spatial variations in spectral properties are averaged out while , as demonstrated by hughes et al . we hypothesized that a more comprehensive modeling of 3d spatial and time dependent transport of bruise chromophores may result in a more accurate determination of the age of the bruise . therefore , our next generation model includes the possibility to model the spatial dependence of the chromophores in the bruise . also , we added horizontal convection and diffusion and michaelis  menten kinetics of the enzymes involved in conversion of hemoglobin to bilirubin . apart from initial formation of the bruise , such a model also needs to account for skin thickness , clearing time of bilirubin , and the different diffusivities of hemoglobin and bilirubin . these different diffusivities imply that the bilirubin ( yellow ) part of bruises diffuses farther than the hemoglobin part . in this article , we present the formulation and parameter analysis of a 3d finite distributed compartment model incorporating all these aspects . we also include a typical natural bruise , photographed every day from 1 to 9  days , and we fitted our model to the measured chromophore areas , as a first quantitative validation of the model . finally , within the validity of our modeling , we show how the age of ( simulated ) bruises can potentially be determined from the time dependent differences between hemoglobin and bilirubin areas within the bruise . the model of the skin consists of three layers ; the top layer of the dermis ( layer 1 ) , the bottom layer of the dermis ( layer 2 ) and the subcutaneous tissue layer ( layer 3 ) . each layer is subdivided into 100    100 corresponding compartments as schematically shown in fig . for ease of analysis , all compartments have the same lateral dimensions equal to the thickness of the total dermis , and the thickness of the layers can be varied to account for inter - individual variability or body location . the initial condition of the bruise is modeled by a pool of hemoglobin of arbitrary shape and size in the subcutaneous tissue ( fig . the bruise develops over time by conversion of hemoglobin into bilirubin and transport of these molecules via pressure driven flow , i.e. convection and concentration driven diffusion . the concentration of hemoglobin within a compartment may change by three processes : ( 1 ) convection in vertical direction from the subcutaneous tissue layer into the dermis ( fig . 1b ) , ( 2 ) vertical diffusion between the layers , as well as horizontal diffusion within the layers ( fig . darcy s law for transport of fluids describes the convection of hemoglobin in vertical direction , assuming no pressure gradient in horizontal direction , hence neglects flow in horizontal direction ( first right hand term in eq .  1 ) . diffusion , in vertical and horizontal direction , follows concentration gradients and is described by the first law of fick ( second right hand term in eq .  1 ) . we will first present the formula for compartments in layer 2 , from which the formulas for layers 1 and 3 follow.fig . 1a the skin model consists of 3 layers , the top layer of the dermis , the bottom layer of the dermis and the subcutaneous fat layer . each layer consists of 100    100 compartmets ( for ease of presentation a smaller number is shown ) . b and c a pool of hemoglobin is defined in the subcutaneous layer . via michealis both hemoglobin and bilirubin flow inside the layers and between the layers a the skin model consists of 3 layers , the top layer of the dermis , the bottom layer of the dermis and the subcutaneous fat layer . each layer consists of 100    100 compartmets ( for ease of presentation a smaller number is shown ) . b and c a pool of hemoglobin is defined in the subcutaneous layer . via michealis both hemoglobin and bilirubin flow inside the layers and between the layers the rate of change of hemoglobin in the compartments in layer 2 ( z  =  2 ) is given by eq . 1):1\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \begin{gathered } { \frac{{\updelta [ { \text{hb}}]_{i , j , z } } } { \updelta t } } = \left ( { { \frac{{[{\text{hb}}]_{i , j , z + 1 } } } { { v_{z + 1 } } } } \cdot { \frac{{k \cdot \updelta p_{z + 1,z } \cdot a_{z + 1 , } } } { { \updelta x_{z + 1,z } } } } } \right ) + \left ( { d_{\text{hb } } \cdot \left ( { \frac{1}{6}\sum\limits_{k = 1}^{6 } { \left ( { { \frac{{\updelta n_{\text{hb } } \cdot a}}{\updelta x } } } \right)_{ijz , k } } + \frac{1}{4}\sum\limits_{l = 1}^{4 } { \left ( { { \frac{{\updelta n_{\text{hb } } \cdot a}}{\updelta x \cdot \sqrt 2 } } } \right)_{ijz , l } } } \right ) } \right ) \hfill \\ - \left ( { { \frac{{v_{\max } \cdot [ hb]_{i , j , z}^ { * } } } { { k_{m } + [ hb]_{i , j , z}^ { * } } } } \cdot [ ho ] \cdot mw_{\text{hb } } } \right ) \hfill \\ \end{gathered } $ $ \end{document } for the compartments in the subcutaneous layer ( layer 3 , z  =  3 ) , eq . 1 is valid if the first right hand term is negated to account for convection out of the subcutaneous layer . neglecting the convection from the subcutaneous layer to the top layer of the dermis ( layer 3 to layer 1 ) , the rate of change for the top layer of the dermis ( layer 1 ) is equal to eq .  1 if the first right hand term is zero . here , [hb]i , j , z is the change in the amount of hemoglobin in compartment i , j , z ( mg ) , 0    i  <  100 , 0    j  <  100 , 1    z    3 , t is the time step ( h ) , [ hb]i , j , z+1 is the amount of hemoglobin in corresponding compartment i , j in the layer z  +  1(mg ) , vz+1 is the volume of the compartment in layer z  +  1 , k is the hydraulic conductivity ( m / nh ) , pz+1,z is the pressure difference between layer z  +  1 and layer z ( n / m ) , az+1,z is the contact surface between the two compartments of layers z  +  1 and z ( m ) , xz+1,z is the distance from center to center of the two compartments in layers z  +  1 and z ( m ) , dhb is the diffusivity of hemoglobin ( m / h ) . to calculate the contributions by diffusion , the average difference in hemoglobin density between compartment i , j , z and the surrounding compartments k and l is used , where k are compartments directly adjacent to compartment i , j , z , and l are compartments across compartment i , j , z , as clarified in fig .  2 . the first sum specifies the difference between compartment i , j , z , and the 6 compartments k , whereas the second sum specifies the difference between compartment i , j , z , and the 4 compartments l. nhb is the difference in hemoglobin density between compartment i , j , z , and k or l ( mg / m ) , a is the contact surface between compartment i , j , z , and k ( m ) and x is the distance from the center of i , j , z , to the center of k ( m ) . for the compartments l , we assumed an effective contact area between compartment i , j , z , and l of a , and a distance from center to center of x2 . vmax is the maximum speed of reaction per mg heme oxygenase ( mol / h / mgho ) , [ hb]i , j , z , is the concentration of hemoglobin in compartment i , j , z , ( mol / l ) , defined as ; \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ [ { \text{hb}}]^ { * } ( { \upmu } { \text{mol}}/{\text{l } } ) = { \frac{{[{\text{hb}}]({\text{mg}}/{\text{l}})}}{{{\text{mw}}_{\text{hb } } ( { \text{mg}}/{\upmu } { \text{mol } } ) } } } $ $ \end{document } , km is the affinity of the enzyme for the reaction ( m =  mol / l ) , [ ho ] is the amount of heme oxygenase in the compartment ( mg ) ( = concentration ( mg / l )    volume ( l ) ) , mwhb is the molecular weight of hemoglobin ; 65  mg/mol.fig . 2the average of adjacent compartments k and across compartments l are used in the calculations the average of adjacent compartments k and across compartments l are used in the calculations for the rate of change of bilirubin , a similar equation as for hemoglobin applies . bilirubin behavior also depends on diffusion and enzymatic conversion ( first and second term in eq .  2 ) . this clearance is not enzyme controlled and is modeled to follow zero order reaction kinetics ( third right hand term in eq .  2 ) . 2):2\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \begin{gathered } { \frac{{\updelta [ { \text{b}}]_{i , j , z } } } { \updelta t}}\,=\,\left ( { d_{\text{b } } \cdot \left ( { \frac{1}{6}\sum\limits_{k = 1}^{6 } { \left ( { { \frac{{\updelta n_{\text{b } } \cdot a}}{\updelta x } } } \right)_{ijz , k } } \,+\,\frac{1}{4}\sum\limits_{l = 1}^{4 } { \left ( { { \frac{{\updelta n_{\text{b } } \cdot a}}{\updelta x \cdot \sqrt 2 } } } \right)_{ijz , l } } } \right ) } \right ) + \left ( { 4 \cdot { \frac{{v_{\max } \cdot [ { \text{hb}}]_{i , j , z}^ { * } } } { { k_{m } + [ { \text{hb}}]_{i , j , z}^ { * } } } } \cdot [ { \text{ho } } ] \cdot { \text{mw}}_{\text{b } } } \right ) \hfill \\ - \left ( { { \frac{{[{\text{b}}]_{i , j , z } } } { { \tau_{\text{b } } } } } } \right ) \hfill \\ \end{gathered } $ $ \end{document}where [b]i , j , z is the change in amount of bilirubin in compartment i , j , z ( mg ) , db is the diffusivity of bilirubin ( m / h ) , nb is the difference in bilirubin density between compartment i , j , z and k or l ( mg / m ) , the number 4 in the second hand term is to account for the fact that 4  mol of bilirubin is formed per mole of hemoglobin , mwb is the molecular weight of bilirubin ; 0.584  mg/mol , [ b]i , j , z is the amount of bilirubin in compartment i , j , z ( mg ) and b is the clearance time of bilirubin from the skin into the lymphatic system ( h ) . for all simulations , the model calculates the spatial distribution of the amount of hemoglobin and bilirubin as a function of time by a forward difference method in space , using labview 8.6 professional development system . the model can simulate a bruise of any given shape and size . for a time step of 0.1  h , a simulation of 400  h takes 1  min on a standard laptop computer . since the injured vessels in the subcutaneous layer close , the convection goes to zero in a short time period , assumed to be 12  h [ 18 , 30 ] . diffusivity values for bilirubin in the skin are unknown , but since 4  mol of bilirubin is produced per mole of hemoglobin , diffusivity values are assumed to be four times larger than hemoglobin diffusivity . the dermis is separated in two layers for two reasons ( 1 ) to enable vertical flow inside the dermis layer , ( 2 ) to enable a comparison between natural bruises and simulated bruises . since only natural bruises in the top layer of the dermis the thickness of the top layer of the dermis was fixed at 400  m , representing the penetration depth of visible light , and the lateral dimensions are equal to the thickness of the total dermis . other parameters ( table  1 ) were varied within biological ranges to account for variability within one individual and between different individuals and to assess the influence of these parameters on the formation and healing of the bruise.table  1parameters used in the modelbiological rangestandard parameters for simulationssi unitstarting concentration of hemoglobin150  g / l150  g / l150  g / lstarting diameter of blood pool2100  mm10  mm0.01  mhydraulic conductivity , k ( t  =  0)5    10  m / nh5    10 m / nh1.4    10 m / nspressure difference , p2.6    10  n / m2.6    10  n / m2.6    10  n / mdiffusivity hemoglobin , dhb1    101    10  m / h1 10  m / h2.8    10  m / sdiffusivity bilirubin , db4    104    10  m / h4 10  m / h1.1    10  m / saffinity , km0.24  m0.24  m0.24    10 mol / lspeed of conversion , vmax3.4  mol / h / mgho3.4  mol / h / mgho9.4    10 mol / s / mghoconcentration of ho0.110  mg / l5  mg / l5    10  g / lclearance of bilirubin , b50400  h150  h54    10  sdermal thickness , x5002000  m , dependent on body site400  m ( top layer)4    10 m ( top layer)600  m ( bottom layer)6    10 m ( bottom layer ) parameters used in the model for the first simulation , a circular pool of hemoglobin was assumed with fixed standard parameters as provided in table  1 . subsequently , we estimated how variation of the clinically most relevant parameters influences the maximal diameter of the bruise in the top layer of the dermis ( for both hemoglobin and bilirubin ) and the time it takes for the bruise to resolve . we assumed a minimal detectable level of chromophores of 1    10  mg per compartment for both hemoglobin and bilirubin . the time point at which the bruise is resolved was taken as the time at which the bilirubin concentration in all compartments is below its assumed critical detectable value . the potential of the model is demonstrated by simulating a natural bruise , originating from a pinch and located on the upper arm of a 26  year old woman ( with assumed skin thickness of 1000  m ) , which was photographed almost every day for 2  weeks . the shape of the non homogeneous starting hemoglobin pool in the subcutaneous layer was estimated from the shape of the bruise in the dermis at day 1 . this shape was obtained using the rgb values of the image ; a small ( hemoglobin containing ) area in the bruise was selected and the range of rgb values in this area was determined . pixels in the image of the total bruise falling within the rgb range of this small area were defined as containing hemoglobin and were given the starting hemoglobin concentration , whereas the other pixels were given a hemoglobin concentration zero . using this method , not only the shape but also the total areas of the hemoglobin and bilirubin parts of the bruise were determined ; the hemoglobin area was obtained by adding up all the hemoglobin containing pixels and multiplying with the area of one pixel . to determine the bilirubin area , a small yellow area was selected . since bilirubin is only produced when hemoglobin is present , the bilirubin area of the real bruise is defined as the hemoglobin area  +  the yellow area . these areas of the real bruise were determined to enable a comparison with the areas of the simulated bruise , which were defined by adding up all compartments having chromophore levels above the detection threshold and multiplying with the area of one compartment . the bruise was then simulated using the hemoglobin containing pixels from the image at day 1 to define the non symmetric shaped blood pool . diffusivity , relaxation time , and concentration of ho were varied until the simulated bruise at various time points resembled the natural bruise at corresponding time points , both in shape as in total area of the hemoglobin and bilirubin parts . an enhanced false color image of the simulation was constructed , in which a high hemoglobin concentration is depicted as a bright red color , fading into a less intense red for lower concentrations , and a high bilirubin concentration is depicted as a yellow color , fading into less intense yellow for lower concentrations . this false color image allows comparison between the simulated hemoglobin and bilirubin areas and the real bilirubin and bilirubin areas , as well as a comparison between high versus low concentrations of the chromophores . the model of the skin consists of three layers ; the top layer of the dermis ( layer 1 ) , the bottom layer of the dermis ( layer 2 ) and the subcutaneous tissue layer ( layer 3 ) . each layer is subdivided into 100    100 corresponding compartments as schematically shown in fig . for ease of analysis , all compartments have the same lateral dimensions equal to the thickness of the total dermis , and the thickness of the layers can be varied to account for inter - individual variability or body location . the initial condition of the bruise is modeled by a pool of hemoglobin of arbitrary shape and size in the subcutaneous tissue ( fig . the bruise develops over time by conversion of hemoglobin into bilirubin and transport of these molecules via pressure driven flow , i.e. convection and concentration driven diffusion . the concentration of hemoglobin within a compartment may change by three processes : ( 1 ) convection in vertical direction from the subcutaneous tissue layer into the dermis ( fig . 1b ) , ( 2 ) vertical diffusion between the layers , as well as horizontal diffusion within the layers ( fig . darcy s law for transport of fluids describes the convection of hemoglobin in vertical direction , assuming no pressure gradient in horizontal direction , hence neglects flow in horizontal direction ( first right hand term in eq .  1 ) . diffusion , in vertical and horizontal direction , follows concentration gradients and is described by the first law of fick ( second right hand term in eq .  1 ) . we will first present the formula for compartments in layer 2 , from which the formulas for layers 1 and 3 follow.fig . 1a the skin model consists of 3 layers , the top layer of the dermis , the bottom layer of the dermis and the subcutaneous fat layer . each layer consists of 100    100 compartmets ( for ease of presentation a smaller number is shown ) . b and c a pool of hemoglobin is defined in the subcutaneous layer . via michealis both hemoglobin and bilirubin flow inside the layers and between the layers a the skin model consists of 3 layers , the top layer of the dermis , the bottom layer of the dermis and the subcutaneous fat layer . each layer consists of 100    100 compartmets ( for ease of presentation a smaller number is shown ) . b and c a pool of hemoglobin is defined in the subcutaneous layer . via michealis both hemoglobin and bilirubin flow inside the layers and between the layers the rate of change of hemoglobin in the compartments in layer 2 ( z  =  2 ) is given by eq . 1):1\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \begin{gathered } { \frac{{\updelta [ { \text{hb}}]_{i , j , z } } } { \updelta t } } = \left ( { { \frac{{[{\text{hb}}]_{i , j , z + 1 } } } { { v_{z + 1 } } } } \cdot { \frac{{k \cdot \updelta p_{z + 1,z } \cdot a_{z + 1 , } } } { { \updelta x_{z + 1,z } } } } } \right ) + \left ( { d_{\text{hb } } \cdot \left ( { \frac{1}{6}\sum\limits_{k = 1}^{6 } { \left ( { { \frac{{\updelta n_{\text{hb } } \cdot a}}{\updelta x } } } \right)_{ijz , k } } + \frac{1}{4}\sum\limits_{l = 1}^{4 } { \left ( { { \frac{{\updelta n_{\text{hb } } \cdot a}}{\updelta x \cdot \sqrt 2 } } } \right)_{ijz , l } } } \right ) } \right ) \hfill \\ - \left ( { { \frac{{v_{\max } \cdot [ hb]_{i , j , z}^ { * } } } { { k_{m } + [ hb]_{i , j , z}^ { * } } } } \cdot [ ho ] \cdot mw_{\text{hb } } } \right ) \hfill \\ \end{gathered } $ $ \end{document } for the compartments in the subcutaneous layer ( layer 3 , z  =  3 ) , eq . 1 is valid if the first right hand term is negated to account for convection out of the subcutaneous layer . neglecting the convection from the subcutaneous layer to the top layer of the dermis ( layer 3 to layer 1 ) , the rate of change for the top layer of the dermis ( layer 1 ) is equal to eq .  1 if the first right hand term is zero . here , [hb]i , j , z is the change in the amount of hemoglobin in compartment i , j , z ( mg ) , 0    i  <  100 , 0    j  <  100 , 1    z    3 , t is the time step ( h ) , [ hb]i , j , z+1 is the amount of hemoglobin in corresponding compartment i , j in the layer z  +  1(mg ) , vz+1 is the volume of the compartment in layer z  +  1 , k is the hydraulic conductivity ( m / nh ) , pz+1,z is the pressure difference between layer z  +  1 and layer z ( n / m ) , az+1,z is the contact surface between the two compartments of layers z  +  1 and z ( m ) , xz+1,z is the distance from center to center of the two compartments in layers z  +  1 and z ( m ) , dhb is the diffusivity of hemoglobin ( m / h ) . to calculate the contributions by diffusion , the average difference in hemoglobin density between compartment i , j , z and the surrounding compartments k and l is used , where k are compartments directly adjacent to compartment i , j , z , and l are compartments across compartment i , j , z , as clarified in fig .  2 . the first sum specifies the difference between compartment i , j , z , and the 6 compartments k , whereas the second sum specifies the difference between compartment i , j , z , and the 4 compartments l. nhb is the difference in hemoglobin density between compartment i , j , z , and k or l ( mg / m ) , a is the contact surface between compartment i , j , z , and k ( m ) and x is the distance from the center of i , j , z , to the center of k ( m ) . for the compartments l , we assumed an effective contact area between compartment i , j , z , and l of a , and a distance from center to center of x2 . vmax is the maximum speed of reaction per mg heme oxygenase ( mol / h / mgho ) , [ hb]i , j , z , is the concentration of hemoglobin in compartment i , j , z , ( mol / l ) , defined as ; \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ [ { \text{hb}}]^ { * } ( { \upmu } { \text{mol}}/{\text{l } } ) = { \frac{{[{\text{hb}}]({\text{mg}}/{\text{l}})}}{{{\text{mw}}_{\text{hb } } ( { \text{mg}}/{\upmu } { \text{mol } } ) } } } $ $ \end{document } , km is the affinity of the enzyme for the reaction ( m =  mol / l ) , [ ho ] is the amount of heme oxygenase in the compartment ( mg ) ( = concentration ( mg / l )    volume ( l ) ) , mwhb is the molecular weight of hemoglobin ; 65  mg/mol.fig . 2the average of adjacent compartments k and across compartments l are used in the calculations the average of adjacent compartments k and across compartments l are used in the calculations for the rate of change of bilirubin , a similar equation as for hemoglobin applies . bilirubin behavior also depends on diffusion and enzymatic conversion ( first and second term in eq .  2 ) . this clearance is not enzyme controlled and is modeled to follow zero order reaction kinetics ( third right hand term in eq .  2 ) . 2):2\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \begin{gathered } { \frac{{\updelta [ { \text{b}}]_{i , j , z } } } { \updelta t}}\,=\,\left ( { d_{\text{b } } \cdot \left ( { \frac{1}{6}\sum\limits_{k = 1}^{6 } { \left ( { { \frac{{\updelta n_{\text{b } } \cdot a}}{\updelta x } } } \right)_{ijz , k } } \,+\,\frac{1}{4}\sum\limits_{l = 1}^{4 } { \left ( { { \frac{{\updelta n_{\text{b } } \cdot a}}{\updelta x \cdot \sqrt 2 } } } \right)_{ijz , l } } } \right ) } \right ) + \left ( { 4 \cdot { \frac{{v_{\max } \cdot [ { \text{hb}}]_{i , j , z}^ { * } } } { { k_{m } + [ { \text{hb}}]_{i , j , z}^ { * } } } } \cdot [ { \text{ho } } ] \cdot { \text{mw}}_{\text{b } } } \right ) \hfill \\ - \left ( { { \frac{{[{\text{b}}]_{i , j , z } } } { { \tau_{\text{b } } } } } } \right ) \hfill \\ \end{gathered } $ $ \end{document}where [b]i , j , z is the change in amount of bilirubin in compartment i , j , z ( mg ) , db is the diffusivity of bilirubin ( m / h ) , nb is the difference in bilirubin density between compartment i , j , z and k or l ( mg / m ) , the number 4 in the second hand term is to account for the fact that 4  mol of bilirubin is formed per mole of hemoglobin , mwb is the molecular weight of bilirubin ; 0.584  mg/mol , [ b]i , j , z is the amount of bilirubin in compartment i , j , z ( mg ) and b is the clearance time of bilirubin from the skin into the lymphatic system ( h ) . for all simulations , the model calculates the spatial distribution of the amount of hemoglobin and bilirubin as a function of time by a forward difference method in space , using labview 8.6 professional development system . the model can simulate a bruise of any given shape and size . for a time step of 0.1  h , a simulation of 400  h takes 1  min on a standard laptop computer . since the injured vessels in the subcutaneous layer close , the convection goes to zero in a short time period , assumed to be 12  h [ 18 , 30 ] . diffusivity values for bilirubin in the skin are unknown , but since 4  mol of bilirubin is produced per mole of hemoglobin , diffusivity values are assumed to be four times larger than hemoglobin diffusivity . the dermis is separated in two layers for two reasons ( 1 ) to enable vertical flow inside the dermis layer , ( 2 ) to enable a comparison between natural bruises and simulated bruises . since only natural bruises in the top layer of the dermis the thickness of the top layer of the dermis was fixed at 400  m , representing the penetration depth of visible light , and the lateral dimensions are equal to the thickness of the total dermis . other parameters ( table  1 ) were varied within biological ranges to account for variability within one individual and between different individuals and to assess the influence of these parameters on the formation and healing of the bruise.table  1parameters used in the modelbiological rangestandard parameters for simulationssi unitstarting concentration of hemoglobin150  g / l150  g / l150  g / lstarting diameter of blood pool2100  mm10  mm0.01  mhydraulic conductivity , k ( t  =  0)5    10  m / nh5    10 m / nh1.4    10 m / nspressure difference , p2.6    10  n / m2.6    10  n / m2.6    10  n / mdiffusivity hemoglobin , dhb1    101    10  m / h1 10  m / h2.8    10  m / sdiffusivity bilirubin , db4    104    10  m / h4 10  m / h1.1    10  m / saffinity , km0.24  m0.24  m0.24    10 mol / lspeed of conversion , vmax3.4  mol / h / mgho3.4  mol / h / mgho9.4    10 mol / s / mghoconcentration of ho0.110  mg / l5  mg / l5    10  g / lclearance of bilirubin , b50400  h150  h54    10  sdermal thickness , x5002000  m , dependent on body site400  m ( top layer)4    10 m ( top layer)600  m ( bottom layer)6    10 m ( bottom layer ) parameters used in the model for the first simulation , a circular pool of hemoglobin was assumed with fixed standard parameters as provided in table  1 . subsequently , we estimated how variation of the clinically most relevant parameters influences the maximal diameter of the bruise in the top layer of the dermis ( for both hemoglobin and bilirubin ) and the time it takes for the bruise to resolve . we assumed a minimal detectable level of chromophores of 1    10  mg per compartment for both hemoglobin and bilirubin . the time point at which the bruise is resolved was taken as the time at which the bilirubin concentration in all compartments is below its assumed critical detectable value . the potential of the model is demonstrated by simulating a natural bruise , originating from a pinch and located on the upper arm of a 26  year old woman ( with assumed skin thickness of 1000  m ) , which was photographed almost every day for 2  weeks . the shape of the non homogeneous starting hemoglobin pool in the subcutaneous layer was estimated from the shape of the bruise in the dermis at day 1 . this shape was obtained using the rgb values of the image ; a small ( hemoglobin containing ) area in the bruise was selected and the range of rgb values in this area was determined . pixels in the image of the total bruise falling within the rgb range of this small area were defined as containing hemoglobin and were given the starting hemoglobin concentration , whereas the other pixels were given a hemoglobin concentration zero . using this method , not only the shape but also the total areas of the hemoglobin and bilirubin parts of the bruise were determined ; the hemoglobin area was obtained by adding up all the hemoglobin containing pixels and multiplying with the area of one pixel . to determine the bilirubin area , a small yellow area was selected . since bilirubin is only produced when hemoglobin is present , the bilirubin area of the real bruise is defined as the hemoglobin area  +  the yellow area . these areas of the real bruise were determined to enable a comparison with the areas of the simulated bruise , which were defined by adding up all compartments having chromophore levels above the detection threshold and multiplying with the area of one compartment . the bruise was then simulated using the hemoglobin containing pixels from the image at day 1 to define the non symmetric shaped blood pool . diffusivity , relaxation time , and concentration of ho were varied until the simulated bruise at various time points resembled the natural bruise at corresponding time points , both in shape as in total area of the hemoglobin and bilirubin parts . an enhanced false color image of the simulation was constructed , in which a high hemoglobin concentration is depicted as a bright red color , fading into a less intense red for lower concentrations , and a high bilirubin concentration is depicted as a yellow color , fading into less intense yellow for lower concentrations . this false color image allows comparison between the simulated hemoglobin and bilirubin areas and the real bilirubin and bilirubin areas , as well as a comparison between high versus low concentrations of the chromophores . a circular symmetric bruise simulated with the standard parameters shows spatial and temporal differences in hemoglobin and bilirubin concentrations . in the center , the kinetics differ from the edge , which has shifted due to diffusion from 5  mm at t  =  0 to 8  mm from the center ( fig .  3 ) . both the center and the edge of the bruise show a fast increase in hemoglobin concentration followed by a fast decrease , and a slower increase and decrease in bilirubin concentration . please note that the peak concentration of hemoglobin in the center is a factor of 6 higher than the peak concentration of hemoglobin at the edge . besides this maximum concentration difference , their temporal behavior is also different : the peak concentration of hemoglobin in the center is reached earlier than at the edge ( center peak at 7.5  h , vs. edge peak at 12  h ) . the peak in bilirubin concentration in the center is also higher than that at the edge . but in contrast to the hemoglobin concentrations , the bilirubin peak in the center is reached later than at the edge ( center peak at 40  h , vs. edge peak at 34  h ) . furthermore , the ratio of hemoglobin over bilirubin in the center is much higher than at the edge.fig . the concentration of hemoglobin in the center in the graph is divided by 2 , to display all chromophores more clearly in 1 graph . the simulations are done with the standard parameters as given in table  1 kinetics in bruise in center and at 8  mm from center . the concentration of hemoglobin in the center in the graph is divided by 2 , to display all chromophores more clearly in 1 graph . the simulations are done with the standard parameters as given in table  1 figure  4 shows the influence of the diffusivity dhb , which is varied between 1    10 and 8    10  m / h for a circular symmetric bruise of a 10  mm starting diameter . we kept db four times larger than dhb and hence dhb is varied between 4    10 and 3.2    10  m / h . a higher diffusivity value ( representing less dense tissue , i.e. less collagen fibers ) leads to a larger bruise ( fig . 4a ) . also , the concentrations of hemoglobin and bilirubin per compartment are lower for higher diffusivity values ( not shown ) and the bruise resolves faster ( fig . m / h , db remained 4    dhb and was varied between 4    103.2    10 the simulations are done with the standard parameters as given in table  1 , except diffusivity a final diameter of bruise for different diffusivities . m / h , db remained 4    dhb and was varied between 4    103.2    10 the simulations are done with the standard parameters as given in table  1 , except diffusivity figure  5 shows the influence of skin thickness on the maximal diameter of a circular symmetric bruise in the top layer of the dermis and the time to resolve . the top layer of the dermis was kept constant at 0.4  mm , and the bottom layer was varied between 0.6 and 1.6  mm . as depicted in fig . 5a , a larger dermal thickness causes a larger diameter of the hemoglobin and bilirubin areas , and a bruise in a thicker dermis takes longer to resolve ( fig . the simulations are done with the standard parameters as given in table  1 , except dermal thickness a final diameter of bruise for different dermal thicknesses . the simulations are done with the standard parameters as given in table  1 , except dermal thickness figure  6a shows the maximal diameter of a circular symmetric bruise for different starting diameters , ranging from 2 to 30  mm . the bilirubin in the bruise diffuses faster than the hemoglobin , resulting in a larger area containing bilirubin than containing hemoglobin . the maximal diameter of the hemoglobin area in the bruise ranges from 8 to 40  mm and the diameter of the bilirubin area ranges from 26 to 96  mm . the concentrations per compartment of hemoglobin and bilirubin in the larger bruises are higher than in the smaller bruises and a larger bruise takes longer to resolve . the diameter of the hemoglobin area of the bruise over time for different starting diameters is shown in fig . the ratio between these two diameters over time for different starting diameters is shown in fig . 6b and c relates the diameter of the hemoglobin area uniquely with the ratio of the hemoglobin over the bilirubin area diameters ( fig . 6d ) , from which the age of the simulated bruise follows directly ; e.g. if the diameter of the hemoglobin area is 20  mm , and the ratio of the hemoglobin over bilirubin area is 0.25 , there is a unique solution for the age of the bruise ; the bruise must be 2.8  days old.fig . b diameter of hemoglobin area of bruise over time for different starting diameters ( the steps in the hemoglobin kinetics are due to the resolution of the model ) . c ratio of the diameter of the area containing hemoglobin over the diameter of the area containing bilirubin for different starting diameters . the simulations are done with the standard parameters as given in table  1 , except for the starting diameter . . ratios of hemoglobin over bilirubin ( color ) a final diameter of the bruise for different starting diameters . b diameter of hemoglobin area of bruise over time for different starting diameters ( the steps in the hemoglobin kinetics are due to the resolution of the model ) . c ratio of the diameter of the area containing hemoglobin over the diameter of the area containing bilirubin for different starting diameters . the simulations are done with the standard parameters as given in table  1 , except for the starting diameter . 6b , c : diameter of the hemoglobin area ( gray ) . ratios of hemoglobin over bilirubin ( color ) for the non homogeneous shaped starting blood pool ( fig . 7b , day 0 ) , derived from the photos of a natural bruise ( fig . 7a , day 1 ) , a qualitative comparison of the temporal and spatial distributions of hemoglobin and bilirubin concentrations in this natural versus the simulated bruise is presented in fig . small amounts of bilirubin are visible . on day 4 , the amount of hemoglobin is lower than on day 1 , and bilirubin has filled the areas between the different hemoglobin parts . on day 7 , the upper part of the bruise has almost resolved , where the lower part contains mainly bilirubin . figure  8 shows the total areas of the hemoglobin and bilirubin parts over time , both for the simulated bruise and for the natural bruise . the total simulated hemoglobin area matches the total hemoglobin area of the natural bruise at several time points for a single set of parameters . the total bilirubin area of the simulated bruise matches the total bilirubin area of the natural bruise in the first 100  h. after 100  h , the total bilirubin area of the simulated bruise is overestimated compared to the natural bruise.fig . b simulations of a non homogeneous bruise ( starting blood pool extracted from photo day 1 ) dermal thickness  =  1000  m , dhb 1.6    10  m / h , db 6.4    10  m / h , b 100  h , concentration ho  =  10.5  mg / l . b total bilirubin area of the simulated and natural bruise a non homogeneous natural bruise photographed on 3 different days ( cross sect . 40  mm ) . b simulations of a non homogeneous bruise ( starting blood pool extracted from photo day 1 ) dermal thickness  =  1000  m , dhb 1.6    10  m / h , db 6.4    10 m / h , b 100  h , concentration ho  =  10.5  mg / l . ballpoint stripes were placed on the same location for comparison a total area of hemoglobin area of the simulated and natural bruise . to our best knowledge , our model is the first that allows simulating the temporal as well as spatial changes in hemoglobin and bilirubin concentrations as occurring during the lifetime of a bruise . particularly , the spatial dependence of the chromophores , as observed in natural bruises , is an essential addition to previous modeling . we showed that combining the hemoglobin area diameter with the ratio of hemoglobin over bilirubin diameters provides a unique solution for the age of circular symmetric simulated bruises ( fig . we hypothesize that this concept can be extended to age determination of inhomogeneous natural bruises because of the faster diffusion of bilirubin compared to hemoglobin , i.e. it fills the gaps between different hemoglobin areas ( fig . . a quantitative validation of the approach was presented , and showed that it was feasible to simulate a non homogeneous bruise with a single set of parameters ( fig . a qualitative validation of this approach requires further study , particularly the assessment of crucial model parameters to match the simulations with the individual bruise under study . we intend to perform this model optimization by combining simulations with spatially resolved reflectance measurements to determine chromophore concentrations   extracted from the spectra using light propagation models [ 4 , 7 , 36 ] . a bruise is formed when small vessels in the subcutaneous tissue layer ( in the order of 510  m ) rupture due to a mechanical impact . although the dermis also contains small vessels , we assume that the blood comes from the subcutaneous layer , therefore , the contribution of ruptured dermal vessels is neglected . we assume a starting concentration of hemoglobin of 150  g / l for all simulations . l for women and 140175  g / l for men   values of children are slightly lower , as are values for anemic patients ) , these small deviations in starting concentration will hardly influence the peak concentration , as well as the temporal behavior of the hemoglobin and bilirubin concentrations . the pressure difference between the blood vessels and the interstitial fluid causes the blood to leak from the ruptured subcutaneous vessels into the surrounding tissue [ 17 , 18 ] . we assume no flow to the epidermal layer because of the tight basal membrane ; therefore this layer is not incorporated in the model . we assume an instantaneously established pool of blood in the subcutaneous layer , meaning no convection horizontally within this layer . the vertical convection from the subcutaneous layer into the dermis is assumed to be slower than the horizontal convection , and is modeled with darcy s law , which uses the hydraulic conductivity k. the value for k is unknown across the subcutis / cutis border , which consists of retinacula , i.e. thin fibrous septa of collagen i  +  iii . in randeberg et al . , a range of values is given for k in the dermis , horizontally along the collagen i  +  iii and elastin fibers , and discussed is how k , vertically across these structures , is smaller . literature shows a broad range of values of k for different types of tissue , and also that k changes with 2 orders of magnitude between normal and edematous skin . this makes it difficult to give an estimate for k. however , because ruptured vessels close ( assumed to be within 12  h ) , the convection gradually goes to zero and the influence of k on the total simulation is small . k for vertical diffusion across the subcutis / cutis border is set to the smallest known value in the dermis at 5    10  m / nh . because of their small influence on the formation and healing , in our simulations k and p are not varied . the free blood initiates an inflammatory reaction , which causes macrophages to extravasate and take up the erythrocytes , free hemoglobin and free heme . the hemoglobin is broken down non - enzymatically into heme and globin [ 8 , 38 ] . inside the macrophages , heme oxygenase-1 ( ho-1 ) converts heme into biliverdin [ 11 , 23 ] . the first step is the rate limiting step and is assumed to occur significantly slower than the second step . menten kinetics , which takes into account the affinity of the enzyme for the reaction ( km ) and the maximum speed of the conversion ( vmax ) . these parameters are standard for each specific enzyme and were , therefore , not varied in this article . the rate of conversion is limited by the enzyme capacity of ho , which is reached even with low concentrations of hemoglobin ( the maximum speed of the reaction is already reached at concentrations a factor of 10,000 lower than the concentration in blood ) . consequently , the kinetics of the hemoglobin and bilirubin do not change for slightly higher or lower starting concentrations of hemoglobin . two other forms of heme oxygenase ( ho-2 and ho-3 ) , both with different enzyme kinetic values , are unlikely to play a role in the bruise and have , therefore , not been taken into account . first , wounds ( e.g. bruises ) show a higher concentration of ho-1 [ 15 , 27 ] and ho-1 is the only ho form that is inducible by high heme concentrations , as occurring in bruises . second , and in contrast , ho-2 is only inducible by adrenal glucocorticoids , implying it will not be upregulated in the bruise . ho-3 has a too low activity in the human body to play a role in bruises [ 10 , 13 ] . the concentration of ho-1 in these simulations is based on the normal serum value of ho-1 , the range is set to 10100 times this normal value , because of the upregulation of ho-1 in wounds . the bilirubin is released from the macrophages into the extracellular space and flows inside the skin before being drained into the lymphatic system . is in the order of 240  h. we varied this clearance time because we hypothesize there are inter - individual variabilities , e.g. more muscle movement will lead to faster lymphatic drainage . we arbitrarily used clearance times between 50 and 400  h , and found this range to affect the kinetics ; a higher clearance time resulted in a longer time to resolve and a larger bilirubin area . bruise , the peak concentration of hemoglobin in the center is higher and is reached earlier than at the edge ( fig . 3 ) , because the hemoglobin diffuses toward the edge and causes a time dependent gradient in concentration . contrastingly , the peak concentration of bilirubin at the edge is reached earlier than at the center , because the pool of hemoglobin at the edge is depleted earlier . in , who assume hemoglobin diffusion to be only in vertical direction , we assume diffusion in both the vertical and horizontal direction . this important addition to the model is based on the fact that diffusion is concentration driven , so irrespective of direction . in addition , the inclusion of horizontal diffusion is the key to the method of age determination described here . diffusivity values of hemoglobin in skin from literature are widespread , as they are adapted from diffusivities in other tissues , e.g. from the diffusivity of myoglobin in skeletal muscle ; 1.08    10  m / h . to account for this uncertainty and also to account for different locations on the body , the diffusivity of hemoglobin dhb was varied between 1    10  m / h and 1    10  m / h , where a higher diffusivity is related to less dense tissue . diffusivity values for bilirubin in the skin are unknown , but since 4  mol of bilirubin is produced per mole of hemoglobin ( as hemoglobin contains 4 heme groups ) , diffusivity values are expected to be larger ; values were varied between 4    10 and 4    10  m / h . this factor between the dhb and db should be verified in future experiments , but our simulation of the natural bruise , where this factor 4 was used , shows at least that this factor gave promising results . the validation of this factor will have to take into account solubility and possible conjugates of bilirubin in the skin . the extent of the mechanical damage on the retinacula , which influences the diffusivity in vertical direction across the subcutis / cutis border , is unknown . however , based on the previous discussion on convective flow across this structure , we assume that diffusivity in the horizontal direction is equal to diffusivity in vertical direction . the diffusivities of hemoglobin and bilirubin are important factors that influence the kinetics of bruises . a higher diffusivity ( corresponding to less dense tissue such as around the eyes ) results in a larger bruise ( fig . also , the higher diffusivity leads to faster kinetics , as the amount of hemoglobin per compartment is reduced and thus less time is needed to convert all the hemoglobin in that compartment to bilirubin . these results imply that the type of tissue influences the rate of development and clearing of bruises . a higher diffusivity may also occur in the case of a severe beating and resulting extensively damaged muscle and collagen fibers . since bilirubin is a small molecule , it diffuses more rapidly in the skin than hemoglobin ( db  >  dhb ) . our assumption that db is not negligible , contrary to randeberg et al . , is an important addition to the modeling of bruises . the assumed faster flow of bilirubin than hemoglobin causes the gaps between the hemoglobin areas to be filled with bilirubin , as shown for the non homogeneous , natural bruise in fig . this phenomenon likely becomes important in future methodology to reliably determine the age of bruises . if the skin thickness used in the model is underestimated compared to the actual thickness ( as simulated in fig . 5 ) , the simulated bruise will resolve faster , resulting in an incorrect prediction of kinetics . the volume of the subcutaneous compartments was kept constant , but as the dermal thickness was varied , the volume and the concentration of the chromophores in the dermal compartments also varied . for a thicker dermis , the concentration is lower , causing a slower flow to the edge of the bruise , and hence smaller maximal diameters of hemoglobin and bilirubin areas ( fig . ideally , the skin thickness should be measured in every situation with high frequency ultra sound or optical coherence tomography . the possible presence of edema , which accumulates in the skin after a traumatic impact and thickens the skin , has not been taken into account in the model . . first , edema increases the skin thickness and consequently induces a slower bruise development and healing . second , the diffusivity is larger because the density of the collagen fibers is reduced due to the increased fluid concentration , resulting in faster kinetics . therefore , the effects of skin edema on the kinetics are expected to be small . increasing the size of the blood pool , as depicted in fig .  6 , results in a larger maximal diameter . interestingly , whereas the maximal diameter of the hemoglobin component of the bruises increases almost linearly with starting blood pool diameter , the maximal diameter of the bilirubin component increases more rapidly for small starting sizes , before increasing linearly . the latter behavior is caused by the lower concentration of bilirubin for the smaller starting sizes ; all the bilirubin becomes cleared before reaching its maximal possible diameter . 6b are due to the resolution of the model ; the lateral dimensions of the compartments in these simulations were 1    1  mm , so if the diameter of the bruise is 10  mm , 10 elements in the cross section of the bruise are above the detection threshold . if the amount of hemoglobin in a compartment becomes lower than the threshold , a step of 1  mm in the computations is seen . the determination of the age of the bruise at 2.8  days can be done under the idealized condition that all parameters are known . if not all parameters are known , a series of measurements with a known time between each measurement can provide enough information to estimate the parameters , as was shown in fig . 8 , and with these known parameters the age of the bruise can be determined . the resulting fit of the parameters for this natural bruise appeared to be within the previously assumed physiological range of parameter values . the shape of the subcutaneous , starting blood pool for the real bruise was estimated using the shape of the visible bruise at day 1 . using the shape at t  =  0 is not possible , because at that time point the blood pool is located in the subcutaneous layer , and the penetration depth of visible light is only 400  m . at day 1 , the visible part of the bruise had already diffused upward from the subcutaneous layer into the upper 400  m   of the skin . although using this diffused part of the bruise as input could lead to an error in the shape and size , the simulation shows good resemblance to the shape of the bruise at day 1 . a good resemblance at the later time points the non homogeneous simulated bruise shows important features that resemble the natural bruise ( fig . 7 ) : first in the shape of the bruise : the bilirubin fills the gaps between the hemoglobin areas in both the simulated and the natural bruise . second , in the high versus low concentrations of chromophores : the amount of hemoglobin is highest on day 1 in both bruises . on day 7 , the upper part of the bruise had almost resolved in both the simulated and natural bruise , and the lower part shows mainly bilirubin . differences between the simulated and the natural bruise exist mainly in the lower part of the bruise ; the amount of bilirubin on day 7 seems higher in the natural bruise . a possible explanation is that during the formation of the bruise , the concentration of hemoglobin in the lower part became higher than the concentration in the upper part . a lower concentration of hemoglobin in the upper part of the bruise could resolve this . another explanation could be the influence of gravity on the flow of the chromophores in the skin , which is not taken into account . a bruise located on the forehead can sag downwards until the patient has one or two bruised eyes . on most parts of the body , where the skin is more dense then around the eyes , gravity is expected to play a less important role . third , the simulated bruise shows resemblance to the natural bruise in total area : the total hemoglobin area of the simulated bruise matches the total area of the natural bruise ( fig .  8) . since extracting the bilirubin area from the photos of the natural bruise proved more difficult than extracting the hemoglobin areas , the deviation from the simulated bilirubin area is larger than for the hemoglobin area . we showed the simulation and age determination of a single bruise . in real abuse situations , where a child is beaten on several days the model does allow addition of an extra ( new ) pool of blood after several days , but the influence of such a second bruise on the age determination requires further study . our model allows simulating the spatial and temporal behavior of the different chromophores in bruises and was used to assess the influence of e.g. diffusivity and skin thickness on the kinetics . varying the parameters in the model in such a way that properties of individual natural bruises are closely matched allows optimizing our model . the spatio - temporal differences of hemoglobin and bilirubin diffusion may identify methodology for more accurate age determination of bruises than currently possible , which could contribute to an improved diagnosis of child abuse .